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A SURVEY  OF  STEADY-STATE  POROUS  FLOW  FREE  BOUNDARY  PROBLEMS 


Colin  W.  Cryer 

0 . Introduction 

A free  boundary  problem  (FBP;  plural  FBPS)  is  a (steady-state) 
boundary  value  problem  involving  differential  equations  on  domains  parts 
of  whose  boundaries,  the  free  boundaries  (FB;  plural  FBS)  are  unknown  and 
must  be  determined  as  part  of  the  solution.  On  such  FBS  the  boundary 
conditions  needed  for  a fixed  boundary  value  problem  (a  problem  for  which 
the  boundaries  are  known)  are  supplemented  by  an  additional  boundary 
condition. 

FBPS  arise  in  porous  flow  problems  when  the  porous  medium  is 
occupied  by  two  fluids  separated  by  a sharp  interface  (the  FB).  The 
most  common  FBPS  involve  water/air  interfaces,  but  water/water  vapor, 
oil/water,  oil/gas,  and  salt  water/fresh  water  interfaces  are  also  often 
considered. 

The  following  terminology  is  often  used  in  porous  flow  FBPS:  FBS 
are  called  free  surfaces  or  phreatic  lines  or  depression  curves  or  water 
tables  or  lines  of  seepage  or  floating  boundaries  or  unknown  boundaries: 
FBPS  are  called  unconfined  flow  problems:  the  porous  medium  is  called 
an  aquifer. 
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Porous  flow  FBPS  have  important  applications  in  the  fields  of  soil 
mechanics,  irrigation,  ground  water  hydrology  (the  control  and  use  of 
water  resources),  petroleum  engineering,  and  industrial  filtration. 

Basic  texts  on  porous  flow  problems  include  those  of  Muskat  [ 19  37  ] 
and  Bear  [ 1972].  More  specialized  tests  include:  Polubarinova  - Kochina 
[ 1 9 62  ] , Harr[l962]  (mathematical  techniques);  Scheidegger  [ i960,  1963], 
Childs  [ 1969  1,  Kirkham  and  Powers  [ 197  2]  (soil  physics);  Todd  [ 19  59  ] 

(ground  water  hydrology);  Carman  [ 19  56],  R.  E.  Collins  [ 1961  ] (gas  flow); 
Luthin  [19  67  ],  van  Schilfgaarde  [ 1974]  (drainage  engineering);  Musicat  [ 1949], 
Pirson  [ 1 9 58 ] (petroleum  engineering);  Hagan,  Haise,  and  Edminster  [ 1967] 
(irrigation);  Remson,  Hornberger,  and  Molz  [ 1 97 1 ] , Oden,  Zienkiewicz, 
Gallagher,  and  Taylor  [ 1 97 4 ] (numerical  methods).  The  survey  paper 
of  Boersma,  Kirkham,  Norum,  Ziemer,  Guitjens,  Davidson,  and  Luthin  [ 1971  ] 
draws  attention  to  important  recent  developments. 

Basic  journals  include: 

Advances  in  Hydroscience, 

J.  Geophysical  Research, 

J.  Hydraulics  Division,  Amer.  Soc.  Civil  Engineers, 

J.  Irrigation  and  Drainage  Division,  Amer.  Soc.  Civil  Engineers, 

J.  Society  Petroleum  Engineers,  Amer.  Institute  Mechanical  Engineers, 

J.  Soil  Science, 
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Water  Resources  Research. 


Despite  the  considerable  literature  on  the  subject,  there  is  no 
comprehensive  survey  of  porous  flow  FBPS.  In  the  present  survey  we 
describe  and  classify  all  the  porous  flow  FBPS  which  have,  to  our  knowledge, 
been  considered  in  the  literature.  It  is  hoped  that  this  will  reduce 
duplication  of  effort  and  also  draw  attention  to  interesting  problems  which 
have  been  neglected. 

In  porous  flow  FBPS  it  is  assumed  that  the  flow  is  saturated,  that 
is,  that  a particular  portion  of  the  porous  medium  is  either  saturated  with 
a fluid  or  contains  none,  so  that  there  are  sharp  interfaces  between  the 
regions  occupied  by  different  fluids.  In  the  context  of  water/air  problems 
the  assumption  that  the  flow  is  saturated  means  that  a portion  of  the  porous 
medium  is  wet  or  dry  but  never  moist.  Porous  flow  FBPS  are  thus  an 
approximation  to  partially  saturated  flow  problems  in  which  the  amount  of 
fluid  in  a particular  portion  of  the  porous  medium  can  vary  between  none 
and  the  maximum  amount  possible.  When  possible  we  draw  attention  to 
the  solutions  of  partially  saturated  flow  problems  which  correspond  to 
particular  porous  flow  FBPS,  because  by  comparing  such  solutions  one  can 
check  the  validity  of  the  assumption  of  saturated  flow. 

The  present  survey  is  one  of  a series  of  surveys  on  FBPS  which  we 
are  writing.  Other  surveys  of  the  series  of  particular  relevance  to  porous 
flow  FBPS  will  be  those  dealing  with  trial  free  boundary  methods  and  with 


variational  inequalities. 


The  remainder  of  this  introduction  is  organized  as  follows:  in 
section  0. 1 we  give  a simple  example  of  a porous  flow  FBP;  in  section  0.  2 
we  discuss  terminology:  in  section  0.  3 we  briefly  describe  some  numerical 
and  analytical  techniques  which  are  mentioned  in  the  remainder  of  the 
survey. 


0.1.  An  example 


The  problem  of  seepage  through  an  earth  dam  (see  Figure  1)  illustrates 


how  porous  flow  FBPS  arise. 


rz 


In  this  problem  water  Irom  an  upstream  reservoir  (or  head  water)  seeps 
through  a rectangular  earth  dam  to  a lower  reservoir  (or  tail  water).  The 
water  only  flows  through  the  region  ft  and  the  remainder  of  the  dam 
remains  dry.  The  air/water  interface  is  denoted  by  r.  It  is  assumed 
that  the  dam  is  so  long  that  the  flow  is  two-dimensional. 

If  Darcy's  law  holds  (see  section  1.1)  then  the  flow  is  described 
by  a velocity  potential  <j>  which  is  related  to  the  hydraulic  head  h by 

4>  = -Kh  , (1) 

where  the  constant  K is  the  hydraulic  conductivity.  In  view  of  (l)  it  is 

immaterial  whether  the  problem  is  formulated  in  terms  of  $ or  h;  here,  we 

usually  use  h.  The  hydraulic  head  satisfies  Laplace's  equation  (see  section  1.1) 

h + h = 0 , in  ft  . (2) 

xx  yy  ’ 

The  boundary  conditions  are  (see  section  2): 

h = H,  on  AB  (interface  with  water  at  rest)  , 

^ - 0,  on  BC  (impervious  boundary)  , 
an 

h = h^,  on  CD  (interface  with  water  at  rest)  , (3) 

h = y,  on  DE  (interface  with  air)  , 

= 0,  on  EA  (streamline)  . 

9n 

If  r were  known  then  equation  (2)  with  boundary  conditions  (3) 
would  suffice  to  determine  h.  Since  r is  not  known  we  have  a FBP 
and  an  extra  condition  is  required  on  r.  This  extra  condition  is: 


h = y,  on  EA  (interface  with  air)  . 


(4) 


This  FBP  is  discussed  further  in  section  3. 1 . 1 . 1 . 


0.2.  Terminology 

As  illustrated  by  the  example  in  section  0.  1,  a porous  flow  FBP 
has  the  following  structure: 


du  = 

o, 

in 

Bu  = 

0, 

in 

Bd  , 

Cu  = 

0, 

on 

r , 

where  u is  the  solution  and  where  d represents  the  linear  or  nonlinear 

elliptic  ditferential  equation(s)  which  must  be  satisfied  by  u on  the 

open  set  P represents  the  linear  or  nonlinear  boundary  conditions  which 

must  be  satisfied  on  the  boundary  3fl  of  fl.  If  a#!  were  completely 

known  then  the  first  two  equations  of  (l)  would  define  u completely. 

However,  part  of  dd,  namely  the  FB  r,  is  unknown  and  C represents 

the  additional  linear  or  nonlinear  boundary  condition  imposed  on  r. 

In  (l)  it  is  to  be  understood  that  the  open  set  A may  consist  of 

m 

the  union  of  several  domains,  A = U jn  , in  which  case  the  first  equation 

i = l 

of  (l)  is  equivalent  to 

d. u = o in  fi.,  1 < i < m . 
l i — 

The  operators  C7.  mayor  may  not  be  the  same. 

The  domains  #S.  will  be  called  free  domains  (FD;  plural,  FDS). 

When  it  is  necessary  the  number  of  domains  will  be  indicated  by  saying 
that  the  problem  is  an  m-domain  FBP.  When  m = 1 we  set  Cl  E d^t  d s 
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Using  this  terminology,  we  may  describe  the  example  in  Section  0.  1 


as  a 1 -domain  FBP  with  FD  fi  and  FB  T. 

Most  of  the  porous  flow  FBPS  which  have  been  considered  are 
1 -domain  FBPS.  A two-domain  FBP  would  arise,  for  example,  in  a coastal 
region  in  which  porous  flow  of  fresh  water  occurred  in  a region  and 

porous  flow  of  salt  water  occurred  in  a region  the  FB  being  the  salt 

water/fresh  water  interface. 

It  is  usually  neither  possible  nor  necessary  to  classify  the  boundary 
conditions  on  the  FB  into  "conventional  boundary  conditions"  and  "supple- 
mentary boundary  conditions",  and  the  division  of  the  boundary  conditions 
on  r into  those  represented  by  8 and  that  represented  by  C is  an 
arbitrary  decision. 

The  FB  r is  defined  to  be  that  part  of  whose  "shape"  is 
unknown;  the  remainder  of  9j0  will  be  called  the  fixed  boundary.  Points  at 
which  the  FB  intersects  the  fixed  boundary  will  be  called  points  of  detach- 
ment. If  the  location  of  a point  of  detachment  is  prescribed  it  is  said 
to  be  fixed;  otherwise  it  is  said  to  be  free.  We  illustrate  these  definitions 
by  considering  Figure  1 of  section  0.  l.  The  point  A is  known  so  that 
A is  a fixed  point  of  detachment.  The  point  E is  not  known  so  that  E 
is  a free  point  of  detachment.  The  FB  is  EA;  although  the  length  of  DE  is 
not  known  its  shape  is  known  so  that  DE  is  part  of  the  fixed  boundary 
and  not  part  of  the  FB. 
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In  figures,  the  notation  of  Figure  1 will  be  used:  FBS  will  be 
denoted  by  broken  lines,  fixed  points  of  detachment  by  solid  circles, 
and  free  points  of  detachment  by  open  circles,  while  the  FD  will  be 
denoted  by  hatching. 

The  FB  r must  usually  satisfy  certain  restraints  which  we  call 
the  auxiliary  restraints.  The  auxiliary  restraints  arise  in  several  ways: 

(i)  In  some  FBPS  the  boundary  operator  C involves  an  unknown 
constant.  In  such  FBPS  it  is  often  the  case  that  the  position 
or  slope  of  one  endpoint  of  r must  be  prescribed  or 
satisfy  an  additional  condition.  FBPS  of  this  type  arise, 

for  example,  when  certain  porous  flow  FBPS  are  reformulated 
using  the  Baiocchi  transformation  (Baiocchi,  Comincioli, 
Magenes,  and  Pozzi  [ 197 3] ) . 

(ii)  It  is  usually  assumed  that  the  solution  u satisfies  certain 
smoothness  conditions.  For  example,  in  the  example  of 
section  0.  l it  is  assumed  that  the  flow  velocity  at  A is 
finite.  It  is  often  a consequence  of  such  assumptions  that 
the  slope  of  r at  its  endpoints  is  either  uniquely  determined 
or  limited  to  a finite  number  of  values;  see  section  2.2.  1. 

(iii)  r must  of  course  satisfy  "obvious"  restraints  such  as  that 
it  should  not  intersect  itself.  While  obvious  in  the  physical 
co-ordinates  such  restraints  become  less  obvious  in  trans- 


formed co-ordinates. 


(iv)  In  some  cases  the  solution  of  a FBP  is  not  unique  and 

auxiliary  restraints  must  be  imposed  in  order  to  eliminate 
undesirable  solutions.  For  example,  Hamel  [ 1938,  p.  43] 
shows  that  the  problem  of  seepage  from  a canal  has  two 
solutions,  one  of  which  is  physically  unacceptable  because 
the  fluid  pressure  is  unbounded. 

0.3.  Brief  summary  of  some  analytical  and  numerical  methods 

In  the  remainder  of  this  survey,  reference  will  be  made  to  some 
analytical  and  numerical  methods.  Here  we  provide  a brief  description 
of  these  methods.  See  also  Magenes  [ 1972]. 

The  Hodoqraph  Method 

There  is  only  one  general  method  for  obtaining  explicit  solutions 
of  FBPS,  namely  the  hodograph  method  (and  variants  thereof).  The  hodo- 
graph  method  relies  heavily  upon  conformal  mapping  and  is  essentially 
restricted  to  FBPS  for  which: 

(i)  The  governing  differential  equation  for  h e * is  Laplace's 

equation  in  the  plane, 

* +4-  = 0,  in  A . 

^xx  ^yy 

(ii)  The  fixed  boundary,  a#-1",  is  polygonal. 

(iii)  The  boundary  conditions  are  such  that  there  are  two  distinct 
conformal  maps  of  A onto  domains  of  known  shape.  In 


particular,  let  ij;  be  the  complex  conjugate  of  <j>  so  that 


f(z)  = + iy)  + iv|/(x  + iy)  is  an  analytic  function  of 

z = x + iy.  Let  w(z)  = df(z)/dz.  Let  and  £ denote 
the  images  of  the  FD  £ in  the  f and  w planes,  respec- 
tively. If  the  boundary  conditions  satisfied  by  h=<j>  are  such 

that  £,  and  £ are  known,  then  the  fact  that  df/dz  = w 
f w 

together  with  the  fact  that  the  conformal  mapping  of  £ onto 

fi  can  be  obtained  makes  it  possible  to  solve  the  FBP 
w 

analytically. 

The  hodograph  method  has  been  extensively  used  with  great  ingenuity 
to  obtain  analytic  solutions.  Harr  [ 1962],  Polubarinova  - Kochine  [ 1962] 
and  Bear  [ 1972]  describe  some  of  the  extraordinary  complicated  FBPS 
which  have  been  solved. 

Solution  in  auxiliary  planes 

In  many  cases  it  is  convenient  to  consider  a FBP  in  an  auxiliary 
plane.  Thus,  with  the  notation  used  to  describe  the  hodograph  method, 
one  may  consider  the  FBP  in  the  f-plane  (or  d4>-plane)  or  the  w-plane 
(or  hodograph  plane)  and  then  solve  numerically  using  finite  differences 
or  other  methods. 

Variational  inequalities 

An  exciting  recent  discovery  due  to  Baiocchi  [ 1 97 1 ] has  been  that 
many  porous  flow  FBPS  can  be  reformulated  as  variational  inequalities. 

In  the  case  of  the  example  of  section  0.1,  let 

H * 

w(x,y)  = f [ h(x,  t)  - t]dt  , 

y 
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h(x,  y)  = 


where 

h(x,  y),  if  (x,  y)  « 0 , 
y , otherwise 
Let  ft  be  the  rectangle  ABCF  (see  Figure  1 of  section  0.1)  and 
let  g be  defined  on  the  boundary  9ft  of  ft  as  follows: 

g = (H  - y)2/Z,  on  AB  , 

= H2/2  + (h2  - H2)x/2f,  on  BC  , 

= (hd  - y)2/2,  on  CD  , 

= 0,  on  DEFA  . 

It  can  be  shown  that  w solves  the  problem:  Minimize  the  integral 

2 . 2 


J(v)  = 7 //  (v  + v + 2v)dxdy 
c ft  y 


among  the  class  of  functions  v which  are  defined  and  non-negative  on 
ft  and  equal  to  g on  3ft.  This  problem  is  one  of  several  equivalent 
formulations  of  variational  inequalities. 

The  use  of  variational  inequalities  has  led  to  profound  results. 

Further  details  will  be  found  in  the  papers  of  Baiocchi  listed  in  the  bibliography. 
Variational  inequalities  can  be  solved  numerically  by  several  methods  including 
finite  differences  and  finite  elements. 

Trial  free  boundary  methods 

We  repeat  the  basic  FBP: 


C7u  = 

o, 

in 

0 , 

* 

Bu  = 

o, 

on 

90  , 

Cu  = 

0, 

on 

r . 
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One  approach  to  solving  this  FBP  numerically  is  to  generate  a sequence 
of  trial  FBS  and  approximate  solutions  uj^  as  follows: 


say. 


Step  0.  Guess  an  initial  trial  FB,  I ^ 

Step  1 ■ Given  let  be  the  corresponding  domain.  Compute 

an  approximation,  uj^  say,  f°  tbe  solution  u^  ^ of  the  problem 

<7u(k)  = 0,  in  jo(k)  , 

Bu^  = 0,  on  • 

Step  2.  Given  and  uj^  compute  a new  trial  FB  r^k  ^ by 

(k) 

requiring  that  Cu^  should  be  approximately  equal  to  zero  on 
r^k+1^;  i.e.  ''move  the  boundary"  from  r^k^  to  r^k  • 
Following  Birkhoff  [ 1961  ] we  will  call  this  method,  and  variations  thereof, 
a trial  free  boundary  method.  (Previously,  we  have  called  methods  of  this 
type  "move-the-boundary " methods  (Cryer  [ 1970a]),  but  Birkhoff's 
terminology  is  better. ) 

, (k) 

In  Step  2 the  computation  of  the  approximate  solution  u requires 
the  numerical  solution  of  a fixed  boundary  value  problem,  which  can  be 
done  using  any  standard  method.  We  then  speak  of  a trial  free  boundary 
method  using  finite  elements,  a trial  free  boundary  method  using  finite 
differences,  etc. 

Cryer  [ 1970a]  gives  a general  discussion  of  trial  free  boundary 
methods. 
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Time-dependent  methods 

Porous  flow  FBPS  have  often  been  solved  b£  numerically  solving 
a time-dependent  porous  flow  problem  and  letting  the  time  become  very 
large,  so  that  steady-state  is  approached.  We  call  such  methods  time- 
dependent  methods,  and  qualify  them  by  indicating  the  numerical  method 
used  to  solve  the  time-dependent  problem. 
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i . The  Governing  Equations 


Most  of  the  TBPS  considered  in  the  literature  concern  the  flow 
subject  to  Darcy's  law  of  an  incompressible  heavy  fluid  through  an 
isotropic  homogeneous  medium,  and  the  reader  should  assume  that  this 
case  is  being  considered  unless  it  is  explicitly  stated  otherwise.  The 
governing  equations  for  this  special  case  are  given  in  section  1.1. 

There  are  minor  but  extremely  irritating  differences  of  notation 
in  the  literature.  The  following  notation  will  be  used  here: 


x = Xj  and  y = 


co-ordinate  vector 


plane  co-ordinates 


r and  y or  x and  y cylindrical  co-ordinates 
q fluid  velocity 


u and  v 


velocity  components 


pa=0 


density  of  fluid 


fluid  pressure 


atmospheric  pressure 


gravity  acceleration  (along  negative  y-axis) 
gravity  vector 


= p/y  = p/pg 


velocity  potential 

stream  function  (in  saturated  flow) 

suction  (in  unsaturated  flow) 


permeability 


f 


permeability  tensor 


k = (k^)  = pK/py 
K = kpg/p 
K = (K^)  = kpg/p 


P 


h = y 


Ar. 

p(P) 


e 

^ = pg 


hydraulic  conductivity 
hydraulic  conductivity  tensor 
dynamic  viscosity 

hydraulic  head  (or  piezometric  head) ; 
this  is  often  denoted  by  v or  * 

moisture  content 

specific  weight  of  water 


It  is  difficult  to  formulate  the  governing  equations  for  porous 
flow  for  several  reasons: 

(i)  The  flow  of  a fluid  through  a porous  medium,  for  example  the 
flow  of  water  through  soil,  is  an  extremely  complex  matter.  The  ground 
consists  of  layers  of  soil  resting  on  layers  of  rock,  and  is  certainly  not 
homogeneous.  Furthermore,  the  different  layers  are  often  not  isotropic. 
(To  mention  but  one  of  the  many  possible  causes  of  anisotropy,  if  a layer 
of  soil  has  been  formed  by  sedimentation  then  flat  particles  tend  to  be 
oriented  with  their  longest  dimensions  parallel  to  the  layer. ) At  the 
microscopic  level  soil  consists  of  a largo  number  of  soil  particles  of 
different  sizes  and  shapes,  the  interstices  between  the  particles  being 
filled  with  water  and  air.  As  water  flows  through  soil  it  is  subject  to 
electrical  and  chemical  processes  which  are  not  fully  understood.  For 
a detailed  discussion  of  soil  physics  see  Childs  f 1969]  , Boar  [ 197  2]  , 
Kirk  ham  and  Powers  [ 197  2]  . 


(ii)  It  ir.  difficult  to  perform  accurate  experiments  for  porous  flow 
problems.  It  is.  for  example,  a non-trivial  matter  to  obtain  representative 
and  reproducible  soil  samples.  Several  years  ago  wc  were  directly 
involved  in  some  porous  flow  experiments  and  were  surprised  at  the 

bulk  of  the  probe  which  was  used  to  measure  the  pore  water  pressure 
p.  In  the  FBI’S  discussed  in  this  chapter  the  FB  is  shown  as  a smooth 
curve.  In  contrast,  experimentally  observed  FBS  are  often  rather 
irregular.  We  have  seen  this  in  the  laboratory,  and  it  can,  for  example, 
also  be  seen  in  the  figures  in  the  fundamental  work  of  L.  Casagrande 
[ 19  32,  19  34]  on  the  flow  through  dams. 

(iii)  The  theory  of  porous  flow  is  not  supported  by  an  underlying  theory 
in  the  way  that,  say,  the  theory  of  heat  conduction  is  supported  by  the 
theory  of  thermodynamics. 

(iv)  There  have  been  many  attempts  to  derive  the  global  (macroscopic) 
governing  equations  by  making  assumptions  about  the  local  (microscopic) 
equations  and  the  structure  of  the  porous  medium  (Bear  [197  2,  p.  161]  , 
Scheidegger  [ 1963,  p.  641]).  In  the  simplest  model  it  is  assumed  that 
the  porous  medium  consists  of  a large  number  of  circular  capillary  tubes 
and  that  the  fluid  is  a Newtonian  fluid.  Much  more  complicated  models 
have  been  considered;  Scheidegger  [ 1966]  gives  a very  readable  concise 
review  and  Bear  [197  2,  p.  161]  gives  more  details  and  describes  recent 
work.  It  seems  fair  to  say,  however,  that  while  these  models  provide 
insight  none  of  them  is  completely  satisfactory. 
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As  a result  of  the  difficulties  mentioned  above  there  is  agreement 
about  the  governing  equations  for  linear  isotropic  homogeneous  saturated 
flow,  but  there  is  not  universal  agreement  about  the  correct  generalizations 
of  these  equations  to  more  complicated  flows. 

The  great  majority  of  the  porous  flow  FBPS  which  have  been 
considered  in  the  literature  are  for  linear  isotropic  homogeneous  saturated 
flow.  Recently,  however,  some  nonlinear,  anisotropic,  or  inhomogeneous 
FBPS  have  been  solved  (numerically);  we  believe  that  such  FBPS  will  be 
increasingly  studied  and  therefore  discuss  their  governing  equations. 

1.1.  Linear  incompressible  saturated  tlow  in  an  icctropic 
homogeneous  medium 

The  first  basic  equation  is  the  equation  of  continuity 

div  a = 0 . 

The  equation  of  continuity  permits  the  introduction  of  a stream  function 
in  plane  and  axially  symmetric  problems: 


(P) 

u = v 

< 

ti 

■ 

-e- 

(A) 

1 . 

u = - 41  , 
x y 

-3- 

-.1  X 

i 

II 

> 

The  second  basic  equation  expresses  Darcy's  law  that  the  flow 
velocity  is  proportional  to  the  gradient  of  the  hydraulic  head: 

» 
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a = -K  grad  h 


where  the  scalar  constant  K is  the  hydraulic  conductivity  of  the  medium. 
Equivalently, 


K k 

q = - — grad(p  + pgy)  = - ” grad(p+pgy)  , 
pg  P 


where  the  scalar  constant  k = p.K/pg  is  the  permeability. 
It  follows  from  Darcy’s  law  that  the  function 


4>  = -Kh  = - - (p+pgy) 
Y P 


is  a velocity  potential: 


a = grad  4, , 


(P) 

c 

11 

-e- 

x_ 

v = a , 

(A) 

u = 4>  , 
Tx’ 

< 

11 

-e- 

< 

In  the  literature,  the  negative  velocity  potential  (or  specific  discharge 
potential)  is  often  used. 

It  follows  immediately  that  in  plane  and  axially  symmetric  problems: 


(P)  $ +4.  = 0, 

Txx  Tyy 


(A)  A + ~ <|>  + u>  = 0 
yxx  x Yx  Yyy 


(P)  4,  + ^ = 0, 

xx  yy  * 


(A)  - — 4,  + = 0. 

xx  x x yy 


These  equations  form  the  basis  of  most  work  on  porous  flow  FBPS. 
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1-2.  Linear  compressible  anisotropic  saturated  ilow 


There  are  three  basic  equations:  the  equation  of  continuity;  the 
equation  of  state;  and  Darcy's  law. 

The  equation  of  continuity  is 

div(p£)  = 0. 

As  in  fluid  mechanics  the  equation  of  continuity  permits  the  introduction 
of  a stream  function  4<  in  plane  and  axially  symmetric  problems: 

(P)  pu  = i|y  pv  = 

(A)  pu  = — pv  = - — iL  . 

x y’  xTx 

The  equation  of  state  expresses  the  density  p as  a function  of 
the  pressure  p.  To  a very  good  approximation,  water  is  incompressible. 
However,  when  studying  the  flow  of  gases  through  soil,  as  is  necessary 
when  considering  natural  gas  wells,  compressibility  must  be  taken  into 
account.  The  equation  of  state  is  usually  taken  to  be  of  the  form 
(Muskat  [ 19  37]) 

P = P(P)  = PQP  e , 

where  pQ,  m,  and  p are  constants.  Particular  cases  are: 

1 
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m * 0,  P 0, 


Liquids,  incompressible: 

Liquids,  compressible:  m = 0,  P > 0, 

Gases,  isothermal:  m = 1,  P = 0, 

Gases,  adiabatic:  m > 1,  P = 0. 

The  linear  form  of  Darcy's  law  which  is  most  often  used  is 

(Scheidcgger  [ 196  3,  p.  636],  Bear  [ 1972,  p.  124,  136]), 

= -K(x)  grad  h 


where 


h 


= y + 


i f _d£_ 
g J p(p) 


and  K(x)  is  a symmetric  tensor,  the  hydraulic  conductivity  tensor, 


K(x)  = 


K11 

KI2 

K13 

II 

K21 

K22 

K23 

-K31 

K32 

K33- 

In  plane  problems  K must  be  of  the  form 


It  should  be  observed  that  there  is  apparently  no  fundamental 
phys.cn)  law  which  would  require  K(x)  to  he  symmetric.  Fliiggc  [1972, 
p.  101]  observes  that  he  tried  to  find  an  example  of  a porous  flow  system 
with  unsy inmetric  K(x)  but  did  not  succeed.  It  may  be  remarked  that  in 
the  case  of  heat  conduction  the  thermal  conductivity  plays  a similar 
role  to  that  o?  the  hydraulic  conductivity  K in  porous  flow.  It  is  a 
consequence  of  the  fact  that  entropy  cannot  decrease,  that  the  thermal 
conductivity  is  symmetric  and  positive  semi-definite  (Truesdell  and 
Toupin[l960,  p.  709],  Carslaw  and  Jaeger  [ 19 59]). 

It  is  of  interest  that  inhomogeneity  is  sometimes  desirable. 

Kealy  and  Soderberg  [ 1969,  p.  30]  make  the  following  recommendation 
for  mill -tailing  dams:  "Never  construct  homogeneous  dams.  Zoned-type 
construction  must  be  used  to  control  seepage  water  and  consequently 
increase  stability". 

1 . 3.  Nonlinear  saturated  flow 

The  governing  equations  discussed  above  assumed  that  there  was 
a linear  relationship  (Darcy's  law)  between  the  gradient  of  the  piezometric 
head,  grad  h,  and  the  velocity  £.  However,  it  has  been  found  experimentally 
that  the  relationship  between  grad  h and  £ is  nonlinear  when  the  flow  rate 
is  very  low  or  very  high.  The  various  nonlinear  generalizations  of  Darcy's 
law  are  described  below. 
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1.3.1.  Lew  speed  How 


It 


has  been  found  that  there  is  very  little  flow  when  J = Igrad  h| 


is  less  than  a constant  JQ,  the  threshold  gradient.  This  may  be  thought  of 
as  being  due  to  intcrrnolecular  attractive  forces  which  are  negligible  at 
higher  flow  rates. 

Various  modifications  of  Darcy's  Jaw  have  been  suggested.  Many  of 
the  modifications  are  of  the  form 


a = 


if  J < J0 


-KF(J0/n  grad  h,  if  J > JQ, 


where  F(t)  is  a given  function.  Among  the  suggestions  for  F(t)  we 
mention: 

F(t)  = 1-t,  (Irmay:  see  Bear  [ 197 2,  p.  127]), 

F(t)  = (l-t)2  - f 1[  (l-t)3/2  arc  tanth  (l-t)1/2  - In  t - (1-t)]  , 

v J 

(Kovacs  [ 1969] ) . 

For  further  references  see  Bear  [ 197  2,  p.  128]  . 

The  concept  of  a threshold  gradient  introduces  a new  FB  namely 
the  FB  separating  the  region  of  flow  where  J < JQ  from  the  region  of  flow 
where  ] > Jq.  So  far  as  we  are  aware,  this  FBP  has  not  been  considered 
in  the  literature. 
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At  high  flow  rates  it  is  found  that  the  flow  rate  is  less  than  that 
predicted  by  Darcy's  law.  This  can  be  attributed  to  the  onset  of 
turbulence  which  causes  increased  energy  losses. 

The  most  general  generalization  of  Darcy's  law  is  that  of 
Engelund: 

grad  h = -F(  !g.l)3. 

where  F(t)  is  a given  function.  Scheidegger  [I960],  Ahmed  and 
Sunada  [ 1969]  , and  Bear  [ 197  2,  p.  182]  summarize  the  various  choices  for 
F (t)  which  have  been  suggested.  Almost  all  the  proposed  functions  are  of 
the  form 

m 1 

Fit)  = a + bt  (Forchheimer) 

where  a,  b,  and  m are  constants. 

An  important  application  of  the  nonlinear  versions  of  Darcy's  law 
arises  in  connection  with  porous  flow  through  rockfill  dams,  that  is  dams 
built  with  rocks  rather  than  soil.  The  flow  through  rockfill  dams  is 
apparently  more  turbulent  than  the  flow  through  corresponding  earth  dams, 
and  it  is  necessary  to  use  a nonlinear  law.  Fenton  [ 1972a]  has  used 

1 Fit)  = btm-1  , 

with  m = 1,2,  ana  1.86.  Volker[l969]  has  used 
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I'll)  . 31V  I 11.  821  t 


and 

.745 

F(t)  = 8.89  3 t ; 

Volkcr  obtained  the  values  of  the  parameters  by  fitting  to  experimental 
data  using  least-squares,  and  he  observes  that  both  forms  of  F (t ) fitted 
the  data  accurately. 

An  interesting  result  of  Volker’s  work  is  that  the  form  of  the  FB 
does  not  depend  very  much  upon  the  choice  of  F(t),  but  that  the  nonlinear 
equations  give  much  more  accurate  results  for  the  discharge  through  the 
dam,  (Volker  [ 1969,  p.  2108  and  p.  2110]). 

Hamel  [ 1934,  p.  156]  discusses  modifications  to  handle  infinite 
speeds  at  points  of  discontinuity. 

1 . 3.  3.  Slip  ilow 

In  low  pressure  gas  flow  it  is  found  that  the  flow  rates  are  higher 
than  theoretically  predicted.  This  is  called  the  slip  phenomenon  or 
Klinkcnberg  effect. 

Klinkenbcrg  has  suggested  that  (Bear  [ 1972,  p.  128]) 
kg  = k^(l  + bfkjJ/p) 

where  k^  and  k^  are  the  permeabilities  to  gas  and  liquids,  respectively. 

If  the  pressure  p varies  substantially  then  this  implies  that  k varies, 

9 

but  we  do  not  know  whether  this  has  been  taken  into  account  in  the 
literature. 
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1.4.  Partially  .saturate :!  flow 


In  the  discussion  above  of  the  flow  of  water  through  soil  two 
assumptions  were  made:  (1)  the  flow  is  saturated,  that  is,  the  soil  is 
either  wet  or  dry;  <ii)  the  motion  of  the  air  in  the  dry  soil  is  negligible. 

In  certain  circumstances  these  assumptions  are  not  reasonable,  and  it  is 
necessary  to  consider  the  movement  of  water  and  air  through  the  soil. 

In  the  general  case  of  multi-phase  flow  the  porous  medium  contains 
m fluids  which  flow  through  the  medium;  this  is  discussed  in  section  1 . 5. 
An  important  special  case  of  two-phase  flow  is  unsaturated  flow  which  is 
often  considered  in  the  water/air  flow  problems  of  irrigation.  In  unsaturated 
flow  problems  the  motion  of  the  air  is  neglected,  but  allowance  is  made  for 
the  presence  of  air. 

Unsaturated  flow  involves  the  introduction  of  a new  variable  G, 
the  volume  of  water  per  unit  volume  of  soil.  For  steady-state  unsaturated 
flow  of  water  the  governing  equations  are  (Braester,  Dagan,  Neuman  and 
Zaslavsky  [1971,  p.  14]): 

div  g.  = 0 , 

= -K  grad  h, 

h = y + p/pg  = y + p/y  , 

where,  in  distinction  to  the  case  of  saturated  flow,  K is  a function  of 
the  present  and  past  states  of  the  flow. 
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It  has  been  found  experimentally  that  K is  a highly  nonlinear 

function  which  exhibits  hysteresis  and  the  relationships  between  K, 

0 , and  ^ = p/y  have  been  stated  in  several  ways.  (E.  E.  Miller  an- 

Klute  [1967],  Rraester,  Dagan,  Neuman,  and  Zaslavsky  [ 1971]  ) . 

In  the  applications  of  interest  to  us,  K and  G have  been  expressed 

as  functions  of  the  pressure  p.  For  p > p =0  K and  0 have  been 

a 

assumed  to  be  constant.  For  p < 0 expressions  of  the  form 


K = 


A(-p)°  + 1 


9 


(1) 

G 

e = 2 

B(-p)P  + 1 

have  been  used,  where  Kq,  Gq,  A,  and  B,  are  constants. 

In  partially  saturated  flows  there  is  no  water/air  FB:  as  can  be 
seen  from  equation  (l),  the  water  content  is  never  zero.  However,  the 
length  of  seepage  surfaces  is  unknown  because  seepage  occurs  only  when 
the  pressure  p is  positive  (see  section  2.  l). 

A simplification  in  partially  saturated  flow  which  is  sometimes 
made  is  that  the  region  of  unsaturated  flow  is  confined  to  a narrow  zone 
above  the  water  table  (p  = 0)  called  the  capillary  fringe  (Bear  [ 197  2; 
pp.  259,  262,  and  480]).  The  term  "capillary  fringe"  is  used  because  the 
fringe  is  thought  of  as  being  produced  by  capillary  forces  (in  analogy  with 
the  rise  of  a liquid  in  a capillary  tube).  In  saturated/capillary  flow  the 
following  assumptions  are  made: 
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(i» 


There  aie  three  regions  - the  satuiated  region,  the  capillary 
region,  and  the  dry  region. 


(ii)  The  governing  equations  are  the  same  in  the  capillary  region  and 
the  saturated  region.  The  dry/capillary  interface  is  a FB.  The 
additional  condition  on  r is  that 


where  the  constant  pc  is  the  capillary  pressure. 

Bear  [197  2,  p.  259]  states  that  when  a capillary  fringe  is  used, 
the  completely  saturated  region  extends  to  the  capillary /dry  interface. 

It  seems  to  us  that  this  approach  can  lead  to  confusion  because  although 
the  governing  equations  are  the  same  in  the  saturated  region  and  the 
capillary  region  the  boundary  conditions  are  different  (see  section  2.l). 

When  comparing  saturated  flows  with  partially  saturated  flows  or 
saturated/capillary  flows,  it  is  reasonable  to  compare  the  water/air  FB  in 
the  saturated  flow  with  the  water  table,  that  is,  the  curve  p = 0,  in  the 
other  flows. 

It  is  reasonable  to  assume  (although  not  proved)  that  as  the 
constants  A and  B in  equation  (1)  tend  to  <*  the  solutions  of  the 
partially  saturated  flow  problem  will  converge  to  the  solution  of 
the  corresponding  saturated  flow  problem.  Thus,  partially  unsaturated 
flows  provide  a means  for  approximating  saturated  flow  FBPS  by  nonlinear 
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fixed  boundary  value  problems.  This  approach  is  related  to  the  idea 
of  approximating  variational  inequalities  by  nonlinear  fixed  boundary 
value  problems. 

For  further  information  about  unsaturated  flow  the  reader  should 
consult  Braester,  Dagan,  Neuman,  and  Zaslavsky  [1971].  freeze  [1971] 
gives  a general  discussion  of  the  numerical  solution  of  the  three- 
dimensional  partially  saturated  problems.  In  addition  we  draw  attention 
to  some  work  on  the  numerical  solution  of  time -dependent  unsaturated 
flows  which  is  not  quoted  by  Braester  et  al:  Amerman  [ 1969]  , 

Carroll  [ 1969] . 

See  also:  Rubin  [ 1968]  , Verma  and  Brutsaert  [ 1970  ] , Jury  [ 197  3 J . 
1.7.  Multi-phase  porous  flow 

In  multi  - phase  porous  flow  the  porous  medium  contains  m > 2 
fluids.  For  example,  in  oil  reservoirs  the  simultaneous  flow  of  oil,  gas, 
and  water  is  encountered.  In  irrigation  problems,  the  water  available  to 
plants  is  in  the  form  of  water  vapour,  so  that  the  two-phase  water/water- 
vapour  problem  is  of  interest  (see  E.  E.  Miller  and  Klute  [ 1967]). 

For  a general  discussion  ee  Bear  [ 197  2,  chapter  9]  . Tor  numerical 
solutions  see  Startzman  [ 1969). 
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1.6.  Electrokinetic  porous  flow 

I 

The  subject  of  electrokinetic  porous  flow  is  concerned  with  the 
flow  of  a fluid  in  a porous  medium  under  the  influence  of  the  hydrodynamic 
and  electric  fields. 

In  plane  problems  the  governing  equations  are  as  follows: 


2h 

ax 


+ k 

ex 


3E 

9x 


ah 

ay 


+ k 

ey 


ay 


(l) 


J = G — + — ~ , 
x x ax  px  ax  ’ 


J = G — + — ^ , 

y y ay  py  ay  ’ 


av 

+ — * =o 
ay 


f 


ax  ay 


* 0, 


(2) 


(3) 


Here,  and  v^  are  the  velocity  components  of  the  fluid,  and 

are  the  current  densities,  h is  the  hydraulic  head,  E is  the  electric 

intensity,  and  K , k , k , k , G , G , p , p , are  constants. 

’ hx’  hy’  ex  ey  x’  y x’  ry’ 

Equations  (1)  are  a modification  of  Darcy's  law,  equations  (2)  are 
a modification  of  Ohm's  law,  and  equations  (3)  express  the  conservation 
of  mass  and  current. 


-29- 


Tho  above  equations  are  given  by  Lewis  and  Garner  [ 197  2]  who  give 
further  references.  Lewis  and  Garner  solve  a fixed  boundary  electrokinetic 
porous  flow  problem  using  finite  elements.  Lewis  and  Humpheson  [ 197  3] 
solve  a time -dependent  electrokinetic  porous  flow  well  problem;  see 
section  3.  4. 1 . 1 . 

1.7.  Plant  roots 

The  flow  of  water  to  plant  roots  is  of  importance  in  agriculture; 
indeed,  it  is  the  main  purpose  of  irrigation  to  provide  an  adequate  flow  of 
water. 

Two  approaches  have  been  used  to  apply  the  theory  of  porous  flow 
to  water  uptake  by  roots  (E.  E.  Miller  and  Klute  [ 1967,  p.  239]): 

(i ) Root  system  model 

The  plant  roots  are  assumed  to  be  distributed  continuously 
throughout  the  soil  and  the  uptake  of  water  is  treated  by  replacing  the 
equation  of  continuity  by  a modified  equation 

div  a + G = 0, 

where  G is  negative  and  represents  the  uptake  of  water  by  the  roots 
per  unit  volume  of  soil. 
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(ii)  Single -root  model 

The  details  of  the  flow  about  a single  root  are  examined.  It  is 
customary  to  assume  that  the  root  is  cylindrical,  and  that  the  water 
flows  through  the  root  wall  at  a given  rate.  The  boundary  condition 
is  thus 

= -K  grad  h = constant, 

on  the  root  wall  x = a. 

It  appears  that  very  little  work  has  been  done  on  this  interesting 
problem.  The  two  approaches  described  above  are  rather  rudimentary. 

A more  sophisticated  approach  would  require  an  understanding  of  the 
mechanism  of  water  absorption  by  roots,  and  Slatycr  [i960]  is  a good 
starting  point  for  a study  of  the  literature  on  this  subject  (see  also 
Hagan,  Haise,  and  Edminster  [1967,  section  VI]). 
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2.  3cun^arv  :oniitlonp 


2.1.  Boundary  zonaltionr:  Fixed  Loundaries 

The  most,  common  fixed  boundary  conditions  are  given  immediately 
below  and  are  illustrated  in  Tigure  1 for  the  case  of  water  seeping  through 
an  earth  dam.  Unless  otherwise  stated,  saturated  flow  is  considered. 


Figure  1:  Common  boundary  conditions  in  porous  flow 


(i)  Impervious  boundary 

An  impervious  boundary  arises  at  the  surface  of  the  rock  layer  or 
retaining  wall.  No  flow  occurs  across  the  boundary  and  the  boundary 
condition  is: 

3^  • n = 0 . 


If  Darcy's  law  holds  then 


k grad  h • ri  = 0. 
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In  problems  for  which  a velocity  potential  <f.  exists,  we  have 


and  in  problems  for  which  a stream  function  vJj  exists,  we  have 


= constant  . 


These  boundary  conditions  hold  for  saturated,  partly  saturated,  and 
saturated/capillary  flow. 


(ii)  Interface  with  fluid  at  rest 

If  the  fluid  at  rest  is  denoted  by  a subscript  2,  then  since  it  is 
at  rest  the  hydrostatic  head  is  constant, 


h 


2 


-dp__ 

p2<p) 


= constant  . 


(1) 


This  is  an  implicit  equation  which  can  be  solved  for  the  pressure 


p2  = P2(V)  . 


At  the  interface  the  pressures  must  be  equal  so  that  the  boundary 
condition  is 


p = Pj  - P2(y), 

or,  equivalently, 
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dp 

Pj(P)  ’ 


pz(yl 

h * hl  * y + g / 

where  a subscript  1,  or  no  subscript,  denotes  the  fluid  in  motion. 

So  far  as  we  are  aware,  boundary  conditions  of  this  generality 
have  not  been  considered  in  the  literature. 


The  constant  h^  is  usually  determined  from  a given  value  of 
p^(v).  In  particular,  if  the  fluid  2 has  an  interface  with  air  at  a height 
y = H then 


dP 

p2(p) 


= H. 


Three  special  cases  have  been  considered: 

(ii ) (a ) Interface  with  gas 

At  a gas  interface  p^fp)  is  so  small  that  the  condition 
h^  = constant  may  be  approximated  to  high  accuracy  by  the 
condition  p^  = constant. 

The  most  common  gas  interface  is  with  air  for  which  we 

have 


so  that  the  boundary  condition  becomes 

P = 0, 
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or,  equivalently , 


h = y. 

(ii ) (b)  Interface  between  incompressible  liquid:; 
If  both  fluids  are  incompressible  then 


h 


2 


y + p2/p29, 


so  that 


P = p29(h2"y)’ 

or,  equivalently, 

p 2 

h = y + — (h  - y). 

P1  2 

(ii)(c)  Water  interface  with  reservoirs,  wells,  head  waters. 
or  tail  waters 

This  is  a special  case  of  (b)  above  with  P|  = = P so 

that 


p = pg(h2  - y)  » 


or,  equivalently, 


h 
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( i i i ) Seepage  surlaco 

When  a FB  approaches  a fixed  impermeable  boundary  or  a downstream 
tailwater  it  cannot  intersect  them.  Instead,  a seepage  surface  forms 
(see  Figure  1).  On  a seepage  surface  the  fluid  is  in  contact  with  the  air 
so  that,  as  in  (ii)(a)  above,  we  have 

p = 0 , 

or,  equivalently, 

h = y . 

In  partially  saturated  flow  a seepage  surface  occurs  when  the  pressure 
p is  greater  than  atmospheric  pressure,  so  that  the  length  of  the  seepage 
surface  is  unknown  a-priori  (see  Freeze  [ 1971a];. 

(iv)  Exposed  capillary  surface  or  unsaturated  surface 

In  a region  of  unsaturated  flow,  or  in  a capillary  fringe,  ® pressure 
p is  less  than  atmospheric  and  it  is  usually  assumed  that  the  fluid  cannot  flow 
across  an  exposed  surface.  Such  exposed  surfaces  are  thus  usually  treated 
as  being  equivalent  to  impervious  boundaries  and  the  boundary  condition  is 

fl  • n = 0. 

For  equivalent  forms  see  (i)  above. 

(v)  Prescribed  heads 

In  certain  problems,  notably  problems  involving  flow  towards 
wells,  it  is  assumed  that  the  head  h is  a given  function  on  part  of  the 
boundary;  usually  h is  assumed  to  be  constant. 
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( v i ) Prescribed  flow  rale 

Sometimes  the  flow  rate  3 • n is  prescribed  on  part  of  the  boundary. 
An  interesting  example  arises  in  the  study  of  the  absorption  of  water  by 
plant  roots;  Gardner  has  analysed  this  problem  on  the  assumption  that  the  root 
absorbs  moisture  at  a constant  rate  (see  E.  E.  Miller  and  Klute  [1967,  p.  2 38  J ) . 

(vii)  Drains 

The  pressure  p is  usually  prescribed  at  a drain  . In  most 
problems  p = p =0. 

(viii)  Semi-pervious  boundaries 

If  two  porous  media  are  separated  by  a thin  layer  of  semi -porous 
material,  then  it  is  often  assumed  that 

SL  • n = (hj  - h2) 

where  h^  and  h^  denote  the  heads  in  the  two  media.  See  Sear  [197  2, 
p.  216]  and  Hantush  [1964]  . 

The  asymptotic  cenavior  of  the  solution  at  the  intersection  of  * 
different  boun  iaries  is  considered  by  Kravtchenko,  Je  Saint-Marc,  and 
Boreli  [ 19  55],  Polubarinova-Kochina  [ 1962,  p.  276]. 
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In  conclusion  we  observe  that  infinite  flow  velocities  occur  in 


the  following  circumstances  (Muskat  [ 1937,  p.  291]  , Bear  [ 197  2,  p.  289]): 

(a)  At  the  point  of  intersection  of  a seepage  surface  with  a water-table. 

(b)  When  the  interface  with  an  upstream  or  downstream  reservoir  has 
a corner  with  interior  angle  (in  the  flow  domain)  greater  than  rr . 

(c)  When  impervious  boundaries,  seepage  surfaces,  or  capillary  exposed 
surfaces  have  corners. 

So  far  as  we  are  aware,  only  Hamel  [1934,  p.  156]  has  explicitly  considered 
the  errors  caused  by  assuming  the  validity  of  Darcy's  law  at  all  flow 
velocities . The  various  high  speed  flow  formulas  discussed  in  section 
1 . 3.  2 automatically  prevent  infinite  velocities. 

2.2.  Boundary  conditions : FBS 

It  is  possible  to  conceive  of  very  complicated  porous  flow  FBPS. 

For  example,  one  could  consider  the  breaking  of  an  earth  dam  duo  to  excess 
pore  water  pressure,  in  which  case  the  dam  faces  would  be  FBS.  However, 
so  far  as  we  are  aware,  the  porous  flow  FBPS  which  have  been  treated  in 
the  literature  involve  only  fluid  interface  FBS  between  two  fluids  denoted 
by  the  superscript  i for  1=1,  2.  Two  basic  boundary  conditions  are 
imposed,  namely  continuity  of  flow  and  continuity  of  pressure: 

(i)  Continuity  of  Flow 

The  following  cases  arise: 

(a)  No  flow.  The  usual  condition  at  a fluid  interface  is 
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that  there  should  be  no  flow  across  the  interface.  The 
boundary  conditions  are  thus  the  same  as  for  an  imparvious 
boundary  (see  section  2.  l)  namely: 


fl  • n = 0 , for  i = 1,  2. 


If  Darcy's  law  holds  then 


grad  h^  • n = 0,  for  i = 1,  2. 


If  velocity  potentials  exist  then 


4>  = 0 for  i = 1,  2 

n ’ ' 


and  if  stream  functions  iV  exist  then 


4*^  = c(1)  = constant,  for  i = 1,  2. 


(b)  Accretion.  It  is  sometimes  assumed  that  there  is 

accretion  on  the  FB,  and  this  is  expressed  by  the  condition 


XI 

a ■ n = N • n, 


where  N is  a specified  flow.  Usually  JW  is  of  the  form 


N = - Nt 


where  ey  denotes  the  unit  vector  in  the  positive  y 
direction,  N > 0 corresponds  to  positive  accretion 
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(such  as  rainfall  or  irrigation)  and  N < 0 corresponds  to 
negative  accretion  (such  as  evaporation).  (Bear  [ 1972, 
p.  256]). 

The  form  of  (1)  is  somewhat  arbitrary.  One  might  perhaps  assume 
that  in  evaporation  the  rate  of  evaporation  is  proportional  to  the  exposed 
area  so  that  (1)  is  replaced  by 

£ • n = -N  > 0 . 

However,  the  boundary  condition  (1)  is  generally  accepted. 

Pozzi  [ 197 -la,  p.  24]  considers  evaporation  from  a dam  and 
observes  that  the  problem  may  not  be  well- posed  if  the  rate  of 
evaporation  (using  (1))  is  too  large. 

(ii)  Continuity  of  pressure 

The  second  condition  at  a fluid  interface  is  that  the  pressure 
should  be  continuous, 


The  following  cases  arise: 

(a)  Interface  with  a fluid  at  rest 

This  is  discussed  in  detail  in  section  2.  l . The  most 
common  cases  are: 

(a)  Water/air  interface  for  which 
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p(1)  = p „ = p = 0; 

water  a 


or  equivalently, 


(P)  Fresh  water/salt  water  interface  for  which 


p'11  • p 


fw 


= p g(h  - y), 
rsw  sw 


or,  equivalently, 


h = h = y + (h  - y), 

fw  p,  sw  1 

fw 


where  the  salt  water  hydraulic  head  h is  a constant. 

sw 


(y)  Capillary  fringe 

A capillary  fringe  is  an  approximation  to  partially 

saturated  flow  (see  section  1.4).  It  is  assumed 

that  the  alr/water  interface  occurs  not  when  the  pressure 

is  zero  but  when  the  pressure  is  equal  to  -p  where 

c 

p >0  is  a constant,  the  capillary  pressure, 
c 


It  is  sometimes  possible  to  combine  the  above  equations  in  a 
useful  way.  Consider  the  conditions  at  an  air/water  interface  with 
accretion,  namely 
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h - y, 


(3) 


q • n = -Ne  • n , 
a.  _ -y  - » 


(4) 


Bear  [1972,  p.  257]  shows  that  if  Darcy's  law  is  of  the  form 

grad  h 


a = - 


ll  0 


22 


then  (3)  and  (4)  lead  to  the  equation 


K .OUT  + K (OH  . N 

Kil  Ux  K22^y 


2K 


22 


\2 

(k  - n\ 

22 

/ 

1 22  / 

K22  • 


Mauersberger  [1965a]  shows  that  if  N = 0 and  Daicy's  law  is  of 


the  form 


q = -K  grad  h 


then  (3)  and  (4)  lead  to 

grad  h • K grad  h - ey  • K grad  h = 0,  (6) 

and  also  Mauersberger  asserts  that  (3)  and  (4)  are  equivalent  to  (3)  and 
(6).  We  believe  that  this  assertion  of  Mauersberger  can  be  rigorously 
Justified  with  the  use  of  the  maximum  principle. 
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2.2.1.  Foinu  of  detachment 

The  behavior  of  a FB  near  a point  of  detachment  is  discussed  by 
Muskat  [ 19  37  , p.  287],  A.  Casagrande  l 1940,  p.  302),  Harr  ( 1962,  p.  2l] 
and  Bear  [ 197  2,  p.  284]  . 

The  analysis  of  the  point  of  detachment  was  first  systematically 
investigated  by  A.  Casagrande  [1940].  There  are  a number  of  subtle  difficulties, 
and  in  the  literature  the  results  of  Casagrande  are  often  quoted  without 
proof,  reference  being  made  to  Casagrande's  paper.  This  is  not  entirely 
satisfactory  since  Casagrande  used  an  engineering  approach.  The  best 
treatment  in  the  literature  is  that  of  Bear  [ 197  2,  p.  284]  who  makes  systematic 
use  of  the  hodograph  plane. 

In  the  Figures  below  we  show  graphically  various  possible  con- 
figurations for  linear  homogeneous  isotropic  saturated  flow;  we  believe 
that  these  are  correct  but,  (as  explained  at  the  end  of  this  section) 

we  are  not  entirely  satisfied  with  all  the  justifications  which  have  been 
advanced  in  the  literature.  For  clarity,  the  curvature  of  the  FB  is 
neglected  in  the  Figures.  The  flow  velocity  and  slope  of  the  FB  at 
the  point  of  separation  are  denoted  by  q and  |3,  respectively.  The 
slope  (3  can  usually  be  found  from  the  condition  that  the  flow  velocity 
must  be  finite  at  the  point  of  detachment. 
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Figure  3: 

Intersection  of  a FB  with 

(a)  an  upstream  wall 

(b)  a downstream  water  table  (no  accretion 


C-><-  <x 


Figure  4:  Possible  effects  with  accretion  (a)  intersection  with  downstream 


water-table,  (b)  intersection  with  downstream  wall 


The  above  Figures  do  not  exhaust  all  the  possibilities.  If  accretion 
is  present  (negative  or  positive)  then  slight  modifications  of  f’iguies  1,  2, 
and  3 are  necessary;  see  Bear  [197  2,  p.  285  and  p.  289]  . 

The  asymptotic  behaviour  of  the  FB  near  a point  of  detachment  ... 
discussed  by  Kravtchenko,  de  Saint-Marc,  and  Boreli  [ 19  5 5],  J.  M.  Taylor 
[ 1 97  1 ] and  Attchison  [ 1972]. 

We  now  consider  some  of  the  shortcomings  of  the  arguments  associated 
with  the  geometries  shown  in  Figures  1 to  4: 

(i)  One  argument  used  to  justify  the  existence  of  seepage  surfaces  is 
that  if  a seepage  surface  did  not  exist  then  the  FB  would  intersect 
the  water-table  and  infinite  velocities  would  occur  (Muskat  [ 19  37  , 
p.  289]  , Bear  [ 1972,  p.  260]).  However,  the  introduction  of  a 
seepage  surface  does  not  prevent  infinite  velocities,  because  the 
flow  velocity  is  infinite  at  the  intersection  of  a seepage  surface  with 
a water  table  (Bear  [197  2,  p.  287],  Muskat  [ 19  37,  p.  291]). 

(ii)  It  is  assumed  that,  in  the  absence  of  accretion,  a FB  slopes  down- 
wards in  the  direction  of  flow.  That  is,  if  the  flow  is  in  the 
positive  x-direction  and  the  FB  is  the  curve  y = y(x)  then 

dy/dx  < 0.  This  is  of  course  a reasonable  assumption  but  the  proof 
(using  the  maximum  principle)  is  seldom  given  (see  Shaw  and 
Southwell  [ 1941,  p.  15]).  This  assumption  is  used  to  exclude 
certain  geometries: 
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(a)  To  explain  why  there  is  a downstream  seepage  surface 

but  not  an  upstream  seepage  surface  (Bear  [ 197  2,  p.  260]). 

(b)  To  exclude  the  possibility  that  at  an  overhanging  Interface 
with  an  upstream  water-table  (figure  la)  the  FB  is  normal 
to  the  interface,  that  is,  p - a ~ v/z. 

(iii)  The  arguments  for  the  existence  of  a vertical  seepage  surface  at 

an  overhanging  interface  with  a downstream  water-table  (Figure  2c) 
are  rather  involved.  The  arguments  of  Bear  f 197  2,  p.  260]  do  not 
hold  in  this  case.  Muskat  [ 1937,  p.  290]  states  that  he  cannot 
justify  the  existence  of  the  seepage  surface  mathematically  except 
in  the  case  when  there  is  no  downstream  water. 

In  conclusion  we  observe  that  we  believe  that  the  behavior  shown 
in  the  above  Figures  is  correct.  However  a clear  statement  of  the  mathematical 
assumptions  involved  is  desirable.  Our  criticisms  are  not  just  a demand 
for  "proofs  of  obvious  facts".  When  FBPS  are  solved  numerically,  it  is 
necessary  to  ensure  that  extraneous  solutions  are  excluded.  The  trial- 
free -boundary  method  has  been  found  to  be  a remarkably 
stable  method  for  solving  porous  flow  FBPS.  We  are  aware  of  only  one 
instance  when  the  trial -free -boundary  method  was  found  to  be  unstable  and 
this  occurred  when  B.  L.  Taylor  and  3rown  [ 1967]  attacker  a problem 
involving  an  overhanging  interface  with  a downstream  reservoir.  We  believe 
that  the  instability  reported  by  Taylor  and  Brown  is  directly  linked  to  the 
mathematical  difficulties  mentioned  in  remark  (iii)  above. 
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3.  Seepage 


There  are  two  basic  types  of  seepage  problem:  seepage  from 
sources  (reservoirs,  ponds,  channels,  canals,  infiltration  ditches); 
seepage  towards  sinks  (wells,  drainage  ditches,  downstream  reservoirs, 
tail  waters).  These  problems  can  be  combined  in  various  ways,  and 
the  geometry  can  be  further  complicated  by  the  introduction  of  impervious 
barriers,  etc.  We  classify  the  problems  as  follows:  flow  through  dams, 
flow  from  channels,  flow  to  drainage  ditches,  flow  towards  wells.  This 
classification  is  purely  one  of  convenience,  and  several  authors  have 
developed  computer  programs  which  can  handle  more  than  one  type  of 
problem. 

R.  L.  Taylor  [1966]  gives  a general  program  for  solving  plane 
and  axisymmetric  porous  flow  problems  with  general  geometry.  Two 
applications  of  this  program  (flow  through  an  inhomogeneous  trapezoidal 
dam,  and  flow  to  a well)  are  given  by  R.  L.  Taylor  and  Brown  [ 1967]  . A 
later  version  of  this  program  is  given  by  Kealy  and  Busch  [1971]  . 

Neuman  [ 197  2]  gives  a general  program  for  solving  plane  and 
axisymmetric  saturated-unsaturated  porous  flow  problems. 

3.  1 . Seepage  througn  jams 

In  many  parts  of  the  world,  Jams  are  constructed  from  earth  or, 
less  frequently,  rock.  Figures  1 and  2 illustrate  some,  but  not  all,  of 
the  various  possibilities. 
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The  material  in  this  section  is  subdivided  as  follows:  simple 


rectangular  dams,  simple  trapezoidal  dams,  polygonal  dams,  and 
other  geometries. 


.ie  through  simple  rectangular  dams 


The  geometry  is  shewn  in  Figure  1 


The  case  of  a homogeneous  dam  with  constant  pressure  (atmospheric-)  on 
the  FB  has  been  considered  by  many  authors  and  serves  as  a model 

problem  for  seepage  through  dams.  This  section  is  therefore  divided 

\ 

into  two  subsections:  the  first  subsection  considers  the  model  problem; 
and  the  second  subsection  considers  generalizations. 


i.  1 . 1 . 1 . Seepage  througn  simple  rectangular  dams:  the  model  proble. 


k 

i 

The  geometry  is  shown  in  Tigure  1 of  section  3.  l . I . 

In  the  model  problem,  which  is  considered  here,  it  is  assumed  that  the 

dam  is  homogeneous  and  that  the  pressure  on  the  KB  is  constant 

(atmospheric). 

The  model  problem  can  be  solved  analytically  using  the  hodograph 
method,  and  there  are  three  distinct  derivations  in  the  literature. 

The  first  analytical  solution  of  the  model  problem  was  given 
by  Davison  in  1932,  and  was  repeated  by  Davison  [ 19 36 , p.  897]. 

The  second  analytical  solution  was  given  by  Hamel  [19  34] ; this 
solution  is  described  by  Muskat  [ 19  37,  p.  303]  and  Bear  [197  2,  p.  398], 
As  is  usual  in  the  hodograph  method,  Hamel's  solution  contains  a number 
of  parameters.  Hamel  carefully  analyzed  his  solution  and  proved 
(Hamel  [19  34,  p.  147])  that  the  parameters  could  be  chosen  so  that  h/H 
assumes  any  value  between  0 and  1 and  i/ H assumes  any  value 
less  than  °°.  In  other  words,  Hamel  proved  that  a solution  exists  for 
all  H > 0,  h < h,  t > 0. 

Hamel  and  Gunther  [ 1935]  and  Muskat  [ 1935]  have  evaluated 
Hamel's  solution  for  several  values  of  the  parameters.  This  work  is 
discussed  by  Muskat  [1937,  p.  309]  and  the  results  are  quoted  in 
Table  1 below. 
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The  tlii r cl  analytical  solution  of  the  model  problem  is  ylven  by 
Polubari  in  iv  i Koch  ire.  | 1962,  p.  28-lJ.  Thi3  solution  gives  the  geometrical 


dimensions  (Polubarinova-Kochina  [ 1962,  p.  289]), 

„ r' ‘ 2 Kf  cri- (n-o)sin2  *1 
l = C J d*  , (l) 

0 vl-a-fp-olsln  'I' 

„ ~ r K[p4(l-p)sin2*]d* 


H = c f NPHi-PJsin  -i'| a* 

0 n/p-oI  (l-p)sin^ 

/ 7 7 

. ^ r~  c Kfa-  sin  'IMsin  ^ d^ 


. „ — r sin  wisin  M'  d'l' 

h = C\  a J /-"2  , 

0 v(l-arsin  ) ( p - ct-  sin  S') 


in  terms  of  parameters  a,  P,  C.  Here,  K denotes  the  complete 
elliptic  integral  of  the  first  kind, 


* r—  r u<P 

K(£)  = K(-v^)  = J -^Y= 

0 vl-£  sin  <7? 


The  height  of  the  seepage  surface  is  (Polubarinova  - Kochina  [1962,  p.  289]), 

h = h + c j/z  K(cos2^-s-ij * ■c-os  * d* (5) 

0 J (1-flj  sin2  vP ) (1  - sin2 

where  = 1 - at  p^  = 1 - p.  The  FB  is  given  in  parametric  form 
(Polubarinova  - Kochina  [ 1962,  p.  290]), 


^ 2 
r K(sin  ^)sin  ^ d'i' 


_ r Ktsin  ) sin  v 

= 1 - C J i 2 2 

0 v(l-asin  ^Kl-Psin  V) 


^ ^ 2 
_ r K(cos  'I'Jsin’I'  d^' 
y = *»+*»  + C / -7-  ^ (7, 

0 v(l-a  sin  *)(1-P  sin  ♦ ) 

for  0 < ^ < n/2.  Polubarinova  - Kochina  [ 1962,  p.  290]  evaluates 


these  formulas  for  various  values  of  the  parameters  «,  p,  C. 


There  a re  two  misprints  In  Polubnrinova-Kochina  [ 1962,  p.  289]  : 

2 

the  term  in  the  denominator  of  (2)  is  given  as  (P-»)  - (1-P)  sin  \P ; the 
range  of  integration  of  (5)  is  given  as  [0,1].  Also,  in  the  integrals  for 
f , H,  and  h given  by  Polubarinova-Kochina  [ 1962,  p.  288]  , the  term 
t,  in  the  denominator  should  be  *J~£. 


In  addition  to  the  three  analytical  solutions  described  above, 
the  model  problem  has  been  solved  approximately  by  many  authors. 


The  results  are  summarized  in  Table 

1 below. 

Author  (Method) 

H 

1 

h 

h 

s 

Hamel  and  Gunther 
[ 193  5]  (Evaluation 
of  Hamel's  solution) 

. 321 ± . 009 

. 1 58  5±  . 002 

.081  ±.00  5 . 

20 5± . 007 

Muskat[  19  3 5]  (see 

. 3223 

. 1619 

.0842 

2067 

Muskat[l9  37,  p.  314]) 

.6695 

.4439 

.1579 

3598 

(Evaluation  of  Hamel's 

.67204 

. 3293 

0 

4300 

solution) 

.8715 

.4843 

0 

5192 

Polubarinova-Kochina 
[1962,  p.  292] 
(Evaluation  of  solution) 

Many  results  plotted 

graphically 

Breitenoder  [ 1942,  p.  56] 
(Graphical  methods  in 
hodograph  plane) 

1. 49 

1.08 

0 

.74 

Wyckoff  a nd  Reed  [ 19  3 5 , 

. 322 

.157 

.08  3 

206  " 

p.  398]  (trial  free 
ary  method;  analog) 

. 68* 

.61* 

0 

29* 

Shaw  and  Southwell 

24 

16 

4 12. 

* 

75 

[1941,  p.  9]  (see  Southwell 
[1946,  p.  207])  (trial 
free  boundary  method; 
finite  differences) 

McNown,  Hsu,  and  Yih  1 2/3  0 .43 

[ 1955,  p.  668]  (trial  free  (24)  (16)  (0)  (10.32) 

boundary  method;  finite 
differences) 


-54- 


Remark  s 


H t h h 

— — — s 


Finnemore  and  Ferry 
[ I960,  p.  1066]  (trial  free 
boundary  method;  finite 
differences)  " 

24 

16 

4 

1 2.  268  ' 

12.  192 
12. 213 

Different 

initial 

guesses 

Herbert  and  Rushton 
[ 1966 ) (trial  free 
boundary  method; 
resistance  network) 

32 

16 

4 

20.8 

] . M.  Taylor  [ 1971  ;p.  56, 

1 

2/3 

0 

. 5268 

(Ax=H/ 24) 

p.  57 , p,  6 5]  and 

(24) 

(16) 

(0) 

(12.64) 

Aitchison  [ 197  2] 

1 

2/3 

1/6 

. 5330 

(Ax-H/24) 

(trial  free  boundary 

(24) 

(16) 

(4) 

(12. 79) 

method;  finite 

1 

2/3 

1/6 

. 5314 

(Ax-H/ 48 ) 

differences) 

(24) 

(16) 

(4) 

(12.75) 

Comincioli,  Guerri, 
and  Volpi  [ 1971] 

. 322 

. 162 

.084 

. 206  1 1 
. 20684 

(Method  M.) 
(Method  M2) 

(trial  free  boundary 
method;  finite 
differences) 

24 

16 

4 

12.892 
1 2. 368 

(Method  Mj) 
(Method  M2) 

Szabo  and  McCaig 
[1968]  (time- 
dependent;  finite 
differences 

19.7 

26.  56 

1.3 

* 

5.  3 

See  France 
et  al.  [ 1971] 

Herbert  [ 1968  ] 
(Time-dependent; 
resistance  network) 

16 

16 

0 

$ 

7.3 

France,  Parekh,  Peters, 

19.7 

26.  56 

1.3 

* 

5.  3 

and  Taylor  [ 1971] 
(time -dependent; 
finite  elements) 

16 

16 

\ 

0 

6.8* 

Comincioli,  Guerri, 
and  Volpi  { 197 1] 
(variational  inequalities; 

. 322 

. 162 

.084 

. 2364 
. 2043 
. 2150 

(Ax=H/l5) 
(Ax-H/ 30) 
(Ax-H/60) 

finite  differences) 

24 

16 

4 

14.  4 

(Ax=H/ 30) 

Kealy  and  Busch  [ 1971, 
p.9l  (trial -free -boundary 
method:  Unite  elements) 
(Also  Kealy  and  Williams 

f 1971,  p.  MS]  ) 

6 

5 

1 

2.6" 

To  hip  1:  Numprlc.il  solutions  for  the  model  problem 

(see  text  U a explanatory  comments) 
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A number  of  explanatory  comments  must  be  made  regarding  Table  1: 


(i)  The  values  marked  by  an  asterisk  have  been  estimated 
by  us  from  graphs. 

(ii)  The  values  in  brackets  have  been  scaled  by  us  from  the 
results  in  the  previous  line  and  rounded  to  four 
significant  figures  for  easy  comparison. 

(iii)  The  results  obtained  using  variational  inequalities  are 

not  strictly  comparable  to  the  other  results  because  a 

different  formulation  is  used.  The  value  of  h shown 

s 

in  Table  1 is  obtained  using  the  first  y-coordinate  for 
which  the  solution  vanishes  on  the  y grid-line  adjacent 
to  the  downstream  face  of  the  dam.  The  results  given  here 
were  obtained  using  a modified  version  of  the  program  of 
Comincioli,  Guerri,  and  Volpi  [ 1971]. 

In  addition  to  the  results  shown  in  Table  1,  both  J.  M.  Taylor  [1971] 
and  Muskat  [ 19  35]  give  the  solutions  for  other  values  of  the  parameters  and 
unpublished  solutions  have  been  obtained  by  Symm  [ 197  5]  (trial  free 
boundary  method  with  integral  equations)  and  Ferriss  [197  2]  (brute  force  trans- 
formation to  a rectangular  domain). 
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S'Atr  * 


J 


The  evaluation  of  the  solutions  of  Hamel  and  Polubarinova- 


Kochina  is  a difficult  task  by  hand  but  relatively  simple  on  a computer. 

It  is  therefore  unfortunate  that  this  work  has  been  overlooked  and  that 
many  authors  have  computed  numerical  solutions  for  the  model  problem 
without  comparing  their  results  with  the  exact  solutions.  In  Table  1 
several  problems  occur  frequently.  For  comparison  we  give  in  Table  2 
the  corresponding  values  obtained  by  us  by  evaluating  formulas  (1)  to  (5). 

The  integrals  were  evaluated  using  a double-precision  version  of  the  CADRE 
subroutine  of  de  Boor  [1971]  . The  values  of  a and  P were  found  using  a 
single-precision  version  of  the  method  of  Brown  [1969]  applied  to  the  equations 


where  t , H,  and  h are  given  by  equations  (1)  to  (3)  and  where  1 , Hf, 

and  h denote  the  desired  values  of  f,  H and  h respectively.  The 
r 

solutions  given  in  Table  2 are  believed  to  be  correct  to  the  number  of 
decimals  shown.  We  are  grateful  to  Julia  Gray  for  her  help  with  the 


programming. 

H l 

h 

h 

s 

a 

£ 

24  16 

4 

12.7132 

. 09  5 1 26 

.465367 

24  16 

0 

1 2.  5674 

0 

.416089 

.322  .162 

.08  4 

. 204460 

. 098  554 

. 1998  23 

Table  2:  Analytical  evaluation  for  the  model  problem 
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Wo  now  consider  the  analysis  of  the  model  problem. 


It  was  found  by  Charni  [ 19  51]  that  the  flow  Q through  the  dam 
is  given  exactly  by  the  formula 

0 = k(H2  - hZ)/2J  , (9) 


where  k denotes  the  permeability.  Formula  (9)  is  derived  by 
Polubarinova  - Kochina  [ 196  2,  p.  282]  and  Bear  [ 197  2,  p.  367], 

Formula  (9)  is  also  obtained  by  using  the  Dupuit  - Forchheimer  approxima- 
tion (see  Bear  [ 1972,  p.  366]).  Many  authors  have  solved  the  model 
problem  numerically,  computed  the  flow  Q,  and  remarked  how  accurate 
the  Dupuit  - Forchheimer  approximation  (9)  is,  not  realizing  that  is  in 
fact  exact. 

There  are  several  proofs  of  existence  and  uniqueness. 

Hamel  [ 19  34]  proved  existence  by  analysing  in  detail  the  analytic 

solution  obtained  using  the  hodograph  method.  Both  Hamel  [19  34, 
p.  147]  and  Davison  [ 1936,  p.  882]  indicated  that  they  believed 

that  the  solution  was  unique,  but  they  were  unable  to  prove  this. 

Miranda  [1 969]  proves  existence  and  uniqueness  using  conformal 

mapping  techniques.  Baiocchi  [ 1971]  showed  that  the  problem 

could  be  reformulated  as  a variational  inequality,  and  proved  existence 

and  uniqueness;  for  a more  detailed  description  of  this  work  see 

Baiocchi  [ 197  2],  Baiocchi,  Comincioli,  Magenes,  and  Pozzi  [197  3]. 
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It  is  also  possible  to  prove  several  results  about  the  FB. 

Davison  [ 19  36a]  used  the  maximum  principle  to  prove  a number  of 
qualitative  properties  of  the  FB.  Shaw  and  Southwell  [ 19 41,  p.  lb] 
use  the  maximum  principle  to  prove  that  the  FB  is  monotonically 
ecreasing.  Miranda  [ 1 9 69 , p.  76],  Baiocchi  [’.972,  p.  124] 

and  Baiocchi,  Comincioli,  Magenes,  and  Pozzi  [ 197  3,  p.  32]  prove 
that  the  FB  is  analytic.  J.  M.  Taylor  [ 1971]  (see  also  Aitchison  [ 197  2] ) 
obtains  an  expansion  for  the  FB  in  the  neighbourhood  of  the  seepage 
face;  this  is  a special  case  of  the  results  ol  Polubarinova-Kochina 
[ 196  2,  p.  276]  and  Kravtchenko,  de  Saint  Marc,  and  Boreli  [ 19  55]  . 

3.  1 . 1 . 2 . Seepage  through  simple  rectangular  Jams:  the  general  case 
The  geometry  is  shown  in  Figure  1 of  section  3.  i . l . 

The  following  generalizations  of  the  model  problem  have  been 
considered: 

Evaporation 

Davison  and  Roscnhead  [1940,  p.  352]  derive  an  analytic 
solution  for  the  case^when  there  is  no  seepage  face.  Pozzi  [ 1974,  p.  14; 
197  4a,  p.  16]  derives  an  analytic  solution  in  three  special  cases,  for  all 
of  which  there  is  no  seepage  face. 

Pozzi  [ 197  4]  reformulates  the  problem  as  a variational  inequality 
and  proves  existence  and  uniqueness  theorems.  Pozzi  [1974a,  p.  6] 
proves  that  the  FB  is  analytic.  Pozzi  [ 197  4a]  also  gives  numerical  results 
obtained  by  solving  the  variational  inequality  using  finite  differences. 
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Infiltration 


Pozzi  [1974a,  p.  24]  obtains  some  non-existence  results  which 
lead  him  to  question  the  correctness  of  the  problem. 

Variable  permeability 

Outmans  [ 1964]  obtains  the  flow  Q for  horizontally  or  vertically 
stratified  dams. 

Baiocchi,  Comincioli,  Magenes,  and  Pozzi  [ 197  3]  prove  existence 
and  uniqueness  for  horizontally  or  vertically  stratified  dams. 

Bend  [ 197  3,  197  4]  uses  variational  inequalities  to  prove  existence 
and  uniqueness  for  the  case 

k = exp(f(x)  + g(y)) . 

3.1.2.  Seepage  through  simple  trapezoidal  dams 

The  geometry  is  shown  in  Figure  1 in  which  a , |3,  H,  h,  and  either 
i or  must  be  specified.  Unless  otherwise  indicated  we  consider 
the  case  of  a homogeneous  isotropic  material. 
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The  case  of  a trapezoidal  dam  is  more  complicated  mathematically 
than  the  case  of  a horizontal  rectangular  dam  because  the  discharge  Q is 
not  known. 

For  the  case  (3  = 90°  Baiocchi,  Comincioli,  Magenes,  and  Pozzi 
[ 1973,  p.  73]  prove  that  there  exists  at  most  one  solution,  and  Comincioli 
[197  4]  proves  that  there  exists  a solution  and  also  describes  a convergent 
numerical  method.  So  far  as  we  know,  the  existence  of  solutions  for 
p \ 90°  is  not  known. 


Numerical  solutions  have  been  obtained  in  the  following  cases: 


Author  (Method) 

_H 

± * 

a. 

6. 

h 

_S 

Wyckoff  and  Reed 
[ 1935,  p.  400] 
(trial  free  boundary 
method;  analog) 

39* 

54* 

* 

1.42  ~ 0 

& 

1.43  * 0 

30° 

45° 

30° 

45° 

* 

. 26 
. 25* 

Yang  [ 1949  , P.  48] 

3 

10  - 0 

45° 

45° 

* 

.75 

Heinrich  and  Desoyer 
[ 1958]  (time-dependent; 
finite  differences) 

8 

44  12  0 

* 

3 

W.D.L.  Finn  [19  67]  (trial 
free  boundary  method; 
finite  elements) 

4 

7.75  - 0 

45° 

60° 

. * 
2.  125 

France,  Parekh, 

12 

42  - 0 

„ -112 
tan  (— 

, -1,12 
) tan  (- 

$ 

) 2.4 

Peters,  and  Taylor 
[1971,  p.  177]  (time- 
dependent;  finite 
elements) 


Fenton  [ 197  2a]  and  Parkin 

1 - .5  .5  33° 

33°  .75 

[1971]  (trial  free  boundary 
method;  finite  elements) 

1 - . 5 • 5 45° 

45°  .75 

Baiocchi,  Comincioli, 
Guerri,  and  Volpi 
[ 1973,  p.  57]  (varia- 
tional inequalities; 

1 3-20 

45°  90C 

finite  differences) 

Comincioli  [ 1974b] 
(variational  and 
quasi  - variational 
inequalities ; 
finite  differences) 

various 

90° 

ft 

Estimated  by  us  from  graphs. 

Table  1.  Numerical  solutions  for  a simple  trapezoidal  dam 


-62- 


Fenton  { 1072a]  lists  his  program  and  gives  results  for  many 
values  of  the  parameters.  The  work  of  ronton  [ 1972a]  and  Parkin  [ 1971] 
is  for  the  general  nonlincai  law  F(t)  - at™  ^ (see  section  1.  3.2), 
but,  for  consistency,  the  results  quoted  above  arc  for  the  linear  case 
m = 1. 

Freeze  [ 1971a]  uses  a time-dependent  method  with  finite 
differences  to  solve  several  problems  of  flow  through  non-homogeneous 
trapezoidal  dams.  Freeze  considers  saturated -unsaturated  flow,  but 
he  compares  his  results  with  those  for  saturated  flow. 

Pettibone  and  Kealy  [1971]  use  the  trial  free  boundary  method 
with  finite  elements  for  flow  through  an  inhomogeneous  trapezoidal  dam 
(a  mill-tailing  dam). 

In  addition,  Weinig  and  Shields  [ 1936]  , A.  Casagrande  [ 1940, 
p.  31 3]  , and  Neuman  and  Witherspoon  [ 1970,  p.  895]  have  obtained 
approximate  solutions,  but  the  numerical  values  of  the  geometrical 
parameters  are  not  given. 

3.1.3.  Seepage  through  polygonal  dams 

There  is  an  immense  number  of  possibilities,  'ihe  following  problems 

have  been  considered: 

Cu t-off  walls 

A cut-off  wall  is  a wall  (usually  sloping  and  impervious)  built  to 
cut  off  the  flow.  Volker  [ 1969]  has  used  the  trial  free  boundary  method 
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with  finite  elements  to  solve  the  problem  of  nonlinear  flow  through  a 
polygonal  dam  with  cut-off  wall  (Figure  1). 


(based  on  Volker  [1969,  p.  2109]) 

Toe  Drains 

Harr  [ 1962,  p.  223]  gives  the  analytic  il  solution  for  a 
trapezoidal  dam  with  toe  drain  (Figure  2). 


Jeppson  [ 1966,  p.  12;  1968]  has  used  finite  differences  in  the 
4>4j-plane  to  solve  the  problem  of  Figure  2 for  several  values  of  H,  L, 
and  a.  feppson  compares  his  results  with  the  graphical  results  of 
A.  Casagrande  [1940]  and  the  analytical  results. 

Shaw  and  Southwell  [1941]  used  the  trial  tree  boundary  method 
with  finite  differences  to  solve  two  problems  of  a trapezoidal  dam  with  toe 
drain  on  a porous  layer  (Figure  5).  In  the  first  problem  the  material  was 

isotropic  and  homogeneous  (Figure  3o);  on  the  second  problem  there  was 
a sublayer  of  different  permeability.  (Tigure  3b).  Neuman  and  Witherspoon 
[ 1970,  p.  896]  use  the  triol-free-boundary  method  with  finite  elements  to 
solve  a similar  problem. 


(a)  Homogeneous  layers 


Figure  3;  Trapezoidal  dams  with  toe  drains  on  pervious 

layers  (based  on  Shaw  and  Southwell  [1941,  p.  4]) 
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Trapezoidal  dams  with  rock  toes 

A trapezoidal  dam  with  rock  toe  is  shown 


in  Figure  4. 


Figure  4.  A trapezoidal  dam  with  a rock  toe 

Harr  [ 1962,  p.  210]  gives  the  analytical  solution  for  the  case 
when  there  is  no  seepage  face. 

Numerical  solutions  for  the  case  when  there  is  a seepage  face 
have  been  obtained  by:  R.  L.  Taylor  and  Brown  [ 1967]  (trial  free  boundary 
method  with  finite  elements);  Neuman  and  Witherspoon  [ 1970]  (trial  free 
boundary  method  with  finite  elements);  and  Neuman  and  Witherspoon  [ 1971a] 
(time-dependent  method  with  finite  elements.  ) 

v. 

Mill-tailing  dams 

Mill-tailing  dams  are  dams  which  are  built  to  contain  industrial 
effluents.  Kealy  and  Soderberg  [ 1969]  give  a general  survey  of  the  design 
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problems  associated  with  mill-tailing  dams.  The  most  important  aspect  is 
the  stability  of  the  dams.  Instability  can  lead  to  disasters  such  as  at  Aberfan 
In  Wales  in  which  a school  was  buried.  The  location  of  the  TB  substantially  affects 
the  stability  characteristics  of  a dam  (Kealy  and  Williams  [1971,  p.  153]. 

Kealy  and  Busch  [ 1971]  (see  also  Kealy  and  Williams  [ 1971,  p.  152]  ) 
use  the  trial  free  boundary  method  with  finite  elements  to  solve  anisotropic 
non-homogeneous  mill-tailing  dam  problems  (Figure  5). 


Figure  5.  A typical  mill-tailing  dam  (based  on 
Kealy  and  Busch  [ 1971,  p.  16]  ). 


An  interesting  feature  of  these  problem  is  that  sometimes  the  FB  is  concave 
(as  in  Figure  5).  This  is  in  agreement  with  observations  (Kealy  and  Busch 
( 1971,  p.  6 and  p.  25). 
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Sheetpiles 


The  problem  of  a dam  with  sheetpiles  introduces  interesting 
additional  mathematical  and  numerical  difficulties  because  the  upstream 
point  of  detachment  is  free.  (See  Tiyure  6). 

Baiocchi,  Comincioli,  Magcnes,  and  Pozzi  [ 1973,  p.  46]  prove 
existence  and  uniqueness  for  a rectangular  dam  (Figure  6).  So  far  as  we 
are  aware,  there  are  no  numerical  solutions. 


Harr  [ 1962,  p.  176]  gives  the  analytic  solution  for  a dam  with 
sheetpile  wall  and  toe  drain  (Figure  7). 

■ 
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sheet  pile 


-4> 


\ 


Figure  7.  Dam  with  sheet  pile  and  toe  drain 

W.  D.  L.  Finn  [ 1967,  p.  46]  uses  the  trial  free  boundary  method 
with  finite  elements  to  solve  the  problem  of  a trapezoidal  dam  with  two 

sheetpiles.  Maione  and  Franzetti  [ 1969]  give  experimental  results  ior 
a particular  dam,  the  Kainji  dam  in  Niger. 

Trapezoidal  dams  with  sloping  bases 


The  geometry  is  shown  in  Figure  8. 


Figure  8.  Trapezoidal  dam  on  a sloping  base 


Baiocchi,  Comincioli,  Magenes  and  Pozzi  [ 1973,  p.  76]  use 
variational  inequalities  to  prove  uniqueness  for  the  case  (3  = v /2,  and 
Baiocchi,  Comincioli,  Guerri,  and  Volpi  [ 1973,  p.  65]  give  numerical 
results  for  the  case  a = p = ff/2. 

For  the  case  p = tt/2,  Comincioli  [ 1974a]  proves  existence 
and  uniqueness  using  a modification  of  the  method  of  variational 
Inequalities.  Comincioli  also  gives  some  numerical  results. 

Harr  [ 1962,  p.  208]  gives  an  exact  analytic  solution  for  the  case 
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General  polygonal  dnms_ 

Davison  [ 1936a]  considers  the  flow  through  general  polygonal 
dams  [ Figure  9). 


Figure  9:  A general  polygonal  dam 

• 

Three-dimensional  dams 

Freeze  [ 1971a]  uses  the  time-dependent  method  with  finite 
differences  to  solve  a problem  of  flow  through  a three-dimensional 
polyhedron. 

Three-dimensional  problems  have  frequently  been  solved  using 
electrolytic  tanks  (Bear  [ 1972,  p.  702]  , Polubarinova-Kochina  [ 1962, 
p.  463]  . 
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The  following  problems  have  boon  considered: 


Herbert  and  Rushton  f 1966,  p.  60]  use  the  trial  free  boundary 
method  with  resistance  networks  for  the  seepage  through  a dam  with 
parabolic  upstream  face  and  drain  (Figure  1).  They  compare  their  results 
with  the  exact  solution  of  Kozeny. 


(based  on  Herbert  and  Rushton  [1966,  p.  68  J) 


Volker  f 1969]  uses  the  trial  free  boundary  method  with  finite 
elements  to  solve  the  problem  of  nonlinear  flow  through  a "trapezoidal" 
dam  with  a slightly  curved  downstream  face. 

Baiocchi  [ 197  4]  has  shown  how  the  problem  of  flow  through  a 
dam  of  general  cross-section  (Figure  2)  can  be  reformulated  as  a quasi- 
variational  inequality.  Comincioli  [ 1974b]  applies  this  method  to 
obtain  numerical  solutions;  Comincioli  gives  numerical  results  for  a 
trapezoidal  dam  but  points  out  that  the  method  can  be  applied  to  a wide 
range  of  problems. 


Figure  2:  A dam  with  general  cross-section 
3.2.  Seepage  from  channels 

There  is  no  standard  terminology,  and  channels  are  also  called 
infiltration  ditches,  canals . streams . or  ponds  (axisymmetric) . 

It  is  worth  mentioning  (Garg  and  Chawla  [ 1970,  p.  1261])  that  the 
analysis  of  seepage  losses  from  channels  is  an  important  problem.  In 
some  irrigation  systems  in  India  so  much  water  is  lost  from  the  channels 
bringing  water  to  the  fields  that  less  than  half  the  water  originally  available 
is  delivered  to  the  fields.  Not  only  is  that  inefficient,  but  the  ground 
alongside  the  channels  becomes  water-logged  which  creates  additional 
difficulties.  Hagan,  Haise,  and  Edminster  [ 1967,  chapter  1 ) give  a very 
interesting  account  of  the  irrigation  systems  of  ancient  civilizations. 

Harr  [1962,  chapter  9]  and  Polubarinova-Kochina  [1962,  chapter  5] 
give  useful  summaries  of  the  available  analytical  solutions. 


One  minor  point  which  does,  however,  deserve  comment  concerns 
the  behavior  of  the  FB  at  a point  of  detachment  in  the  case  of  a capillary 


fringe.  In  the  literature  the  FB  is  often  shown  normal  to  the  channel 
surface  as  at  point  A in  Figure  1 (Polubarinova-Kochina  [ 1962,  p.  160]  , 

Bruch  [ 1969,  p.  1404]).  However,  consideration  of  the  flow  in  the  hodograph 
plane  (Bruch  [ 1969,  p.  1405],  Bear  [ 197  2,  p.  28  5] ) shows  that  the  FB 
is  horizontal  as  at  point  B in  Figure  1. 


Figure  1:  Behavior  of  the  FB  for  a capillary  fringe 


3.2.1.  Single  channei  into  half-plane 

The  geometry  is  shown  in  Figure  1;  water  seeps  from  a channel 
into  the  surrounding  soil. 
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Figure  1 : Seepage  from  a single  channel  into  a half-plane 


The  case  which  has  been  considered  most  frequently  is  the  case 
of  a channel  of  trapezoidal  cross-section  (Figure  2). 


Figure  2:  Seepage  from  a plane  trapezoidal  channel  into  a 

half -pin  no  (with  capillary  fringe) . 


\ 

i 
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Special  cosor.  include: 


(i)  Symmetric  trapezoidal  channel:  a - p 

(ii)  Triangular  channel:  b = 0 

(iii)  Symmetric  triangular  channel:  a = p and  b = 0 

(iv)  Vertical  slit:  a - p = tt/2  and  b = 0, 

Wedernikow  [ 19  37]  obtained  the  analytical  solution  for  a symmetric 
trapezoidal  channel  (including  rectangular,  triangular,  and  slit  channels 
as  special  cases);  see  Harr  [ 1962,  chapter  9]  , Polubarinova-Kochina 
[1962,  page  149].  Later,  Wedernikow  [ 1940]  obtained  the  analytical 
solution  for  a symmetric  trapezoidal  channel  in  the  case  of  a capillary 
fringe;  see  Polubarinova-Kochina  [ 1962,  p.  160]. 

There  are  a number  of  exact  analytical  solutions  for  channels  with 
a curved  "almost  triangular"  cross-section  (Harr  [ 1962,  p.  231]). 

Hamel  [ 19  38]  gives  a series  solution  for  the  case  when  the  channel 
surface  is  of  the  form 

00  00 

x = X a cos(2n+l)w,  y = bw  + Yj  a sin(2n+l)w  . 
n=0  n n=0 

Hamel  [ 19  38,  p.  43]  points  out  that  it  is  necessary  for  uniqueness 

to  assume  that  the  pressure  is  bounded  at  infinity. 

So  far  as  we  are  aware,  there  are:  (i)  no  numerical  solutions; 

(ii)  no  analytic  solutions  for  plane  unsymmetric  trapezoidal  channels; 

(iii)  no  analytic  solutions  for  axisymmetric  problems.  There  are  of  course 
numerical  solutions  for  flow  into  a finite  layer  (see  section  3.2.2). 


r 


S 


The  geometry  is  shown  in  figure  1 for  a trapezoidal  channel: 


water  seeps  from  the  channel  into  the  surrounding  soil  of  finite  depth; 
at  depth  D there  is  a drain  or  pervious  layer. 


We  consider  the  plane  and  axially  symmetric  cases  separately. 

The  case  of  a plane  symmetric  triangular  channel  has  been 
considered  by  several  workers.  Bruch  and  Street  [ 1967]  give  an 
analytical  solution  and  find  that  it  compares  well  with  experimental 
results  for  a - 45°,  H = 4,  D = 46;  Bruch  and  Sainz  [197  2]  use  an 
interactive  computer  terminal  to  evaluate  this  solution.  Jeppson  [ 1968a] 
uses  finite  differences  in  the  44-plane  and  compares  his  results  with 
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those  of  Bruch  and  Street  [1967]  for  the  case  a = 45°,  H = 9,  D = 51; 
Bruch  and  Sainz  [1972]  plot  their  results  (see  Figure  2)  and  compare  them 
with  those  of  Jeppson  (not  shown).  Jeppson  [1968a]  also  gives  results 
for  several  other  values  of  the  parameters. 


Figure  2:  Seepage  from  a triangular  channel  for 

a = 45°.  H - 9.  D = 51  (based  on  Bruch 
and  Sainz  [1972,  p.  522]). 

Plane  symmetric  trapezoidal  channel  problems  have  been  solved 
for  a variety  of  governing  equations  by  Jeppson  using  finite  differences 
in  the  4>i(j-plane: 
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HcCeroncc 
Jeppson  [ 19 68 a ] 
Jeppson  [ 1969 a] 
Jeppson  [ 1968b] 


Jeppson  and  Nelson  [ 1970] 


Governing  Equations 

Homogeneous,  isotropic. 

Isotropic,  k - Ax  + By  + C. 

Two  soil  layers  1 and  2 with  constant 
permeability  tensors. 


If  the  interface  between  the  layers  is  not 
horizontal  then  we  must  have 

k(1)A(1)  = k(?)A<2)  . 

x y x y 
Partially  saturated  flow  with 

kQ/(p/p0)a  + b,  if  p < 0, 

kQ/b,  if  p > 0. 


Jeppson  [ 1968c]  has  considered  axisymmetric  problems.  Jeppson 
uses  finite  differences  in  the  axially  symmetric  4>«4j -plane  and  allows  the 
permeability  to  depend  linearly  upon  the  depth  y.  The  channels  are  not 
specified  but  appear  as  part  of  the  solution.  In  the  problems  solved  the 
channels  are  approximately  trapezoidal,  and  Jeppson  [1968c,  p.  1279] 
asserts  that  the  solution  for  a channel  of  prescribed  shape  could  be 
obtained  by  an  iterative  procedure.  Neuman  and  Witherspoon  11970]  use 
the  trial  free  boundary  method  with  finite  elements  to  solve  the  same  problem 
and  compare  their  results  with  those  of  Jeppson. 
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Bruch  [ 196')]  has  applied  the  hodograph  methoii  to  the  problem 


of  a periodic  array  of  para  del  (plane)  symmetric  triangular  channels 
with  a capillary  fringe  underlain  by  a drainage  layer  at  finite  depth 
(Figure  1).  We  are  not  aware  of  any  numerical  solutions. 


Figure  1:  Seepage  from  periodic  triangular  charnels 

i.  3.  Drainage 

The  objectives  of  agricultural  drainage  are: 

(i)  The  orderly  removal  of  water  that  is  excess  to  economic  crop 
production. 

(ii)  The  establishment  and  maintenance  of  a proper  salt  balance 
in  the  soil,  since  under  certain  conditions,  and  especially  in 
arid  countries,  inadequate  drainage  can  cause  excessive  con- 
centrations of  salt. 
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The  importance  of  drainage  is  indicated  by  the  fact  that  it  is 
believed  that  inadequate  drainage  contributed  significantly  to  the  decline 
of  many  early  civilizations  and  is  a threat  to  many  present  irrigation 
schemes  (Hagan,  Haise,  and  Edminster  [ 1967 , p.  988]). 

Basic  references  on  drainage  include  Luthin  [ 19  57]  , van 
Schilfgaarde  [197  0,  197  4]. 

Drainage  FBPS  consist  of  combinations  of  the  following  components: 

(a)  Sources : Surface  water  applied  either  by  rainfall  or  sprinkler; 
artesian  water  rising  from  underground;  canals;  "foreign  water" 
flowing  from  a catchment  area. 

(b)  Sinks:  Open  ditches;  mole  drains  (unlined);  tile  drains  (lined). 

In  the  case  of  drains  it  is  usual  to  treat  the  drains  as  line  sinks; 

however,  drains  of  finite  cross -section  are  sometimes 
considered.  Tile  drains  ore  lined  with  a porous  material  which 

introduces  an  additional  resistance  to  flow  (usually  neglected). 

Geometry:  The  two  basic  geometries  are  shown  in  Figures  1 and  2. 

In  Figure  1 a scries  of  parallel  drains  is  fed  by  a uniform  rainfall. 

In  the  case  of  line  sinks,  infinite  depth,  with  or  without  a capillary 
fringe  an  analytical  solution  can  be  obtained  by  the  hodograph  method 
(see  Childs  [ 1959],  van  Schilfgaarde  [ 1970,  p.  65]).  The  case 
of  line  sinks  in  soil  of  uniform  depth  (as  in  Figure  1)  has  been  solved 
approximately  using  infinite  series  by  Kirkham  [1966]  , and  List  [1964] 
gives  another  approximate  solution. 
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Childs  has  used  the  trial -free -boundary  method  with  conducting 
paper  to  study  several  parallel  drain  problems  for  soil  of  uniform  depth: 
Childs  [ 1943]  considers  circular  drains  of  finite  diameter:  Childs  [ 1945] 
considers  (a)  the  effect  of  slits  made  to  observe  the  water  table,  (b)  the 
relative  efficacy  of  flooded  and  empty  drains,  (c)  the  effect  of  putting  a 
drain  at  the  bottom  of  a ditch  and  filling  in,  (d)  the  role  of  the  capillary 
fringe;  Childs  [ 1945a]  considers  the  validity  of  the  assumption  that  the 
flow  of  rainwater  is  uniform  across  the  FB.  Van  Deemter  [ 19  50]  solved 
several  parallel  drain  problems  using  the  trial  free  boundary  method  with 
finite  differences. 

Figure  2 shows  a drain  across  a sloping  bed  down  which  "foreign 
water"  flows  at  a steady  rate.  Childs  [1946]  uses  the  trial-free-boundary 
method  with  conducting  paper  to  study  this  problem  for  both  tile  drains 
and  open  ditches. 


i ! I 


4 
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riguro  2:  Foreign  water  draining  down  a slope 
We  conclude  with  some  remarks: 

(i)  Drainage  problems  involving  drainage  into  ditches  are  equivalent 
to  darn  problems  (see  section  3.1).  Problems  involving  ditches 
are  also  considered  in  section  3.  2. 

(ii)  Drainage  FBPS  have  been  relatively  little  studied  in  comparison 
with  problems  such  as  porous  dam  FBPS,  and  there  is  scope  for 
further  work  such  as  the  application  of  variational  inequalities. 

5.  4.  Seepage  towards  wells  or  ditches 

Wells  play  an  important  role  in  the  use  of  water.  The  books  of 
Harr  [ 1962,  chapter  10]  and  Polubarinova-Kochinn  [ 1962,  chapter  9] 
provide  much  information,  and  Hantush  [ 1964]  gives  a lengthy  survey. 

Wells  also  play  an  important  role  in  the  recovery  of  oil  and  gas, 
but  these  problems  involve  two  fluids  and  are  discussed  in  section  4. 
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Seepage  towards  u .ingle  well 


If  the  well  extends  to  the  bottom  of  the  porous  aquifer  it  is  said 


to  be  fully  penetrating;  otherwise,  it  is  said  to  be  partially  penetrating. 
These  two  types  of  wells  are  considered  separately  below. 


3.4. 1.1.  Seepage  towards  a single  fully  penetrating  well 

The  geometry  is  as  shown  in  Tigure  1.  Water  seeps  towards  a 

well  of  radius  r . The  water  is  pumped  out  at  the  steady  rate  so  that  the 
w 

height  of  the  water  in  the  well  remains  h . The  height  of  the  seepage 

surface  is  h . The  catchment  area  of  the  well  is  of  radius  r and  height 
s e 

h . The  quantity  h - h is  called  the  draw-down, 
e e s 

The  boundary  conditions  on  two  of  the  boundaries  require  comment: 

(i)  The  water  drawn  from  the  well  must  be  balanced  by  a flow 

into  the  region  at  the  outer  boundary  r = r . It  is  usual  to  assume 
that  h = hc  at  r = rg;  this  is  equivalent  to  assuming  that  the 
catchment  area  is  surrounded  by  water  so  that  the  well  is  the  axially 
symmetric  equivalent  of  flow  through  a rectangular  dam. 

(ii)  At  the  well  face  r = r it  is  usually  assumed  that  h = h if 

w w 

y < h arid  h = y if  y > h . However,  wells  are  often  lined  by 

a porous  material  which  reduces  the  flow,  and  Hall  [ 1955,  p.  30] 

allows  for  this  by  setting  h = y + h • on  the  seepage  face  where 

max 

h is  a positive  constant, 
max 


The  subsections  below  consider  the  homogeneous  and  inhomogeneous 


cases,  respectively. 


■ -'T-'.l- 

■ YKf  r*  *V. 


single  fully-penetrating  well 


I 


Here  we  consider  a single  fully  penetrating  well  In  a homogeneous 


isotropic  medium  with  permeability  k.  The  geometry  is  shown  in  Figure  1 
of  section  3.4.1.  l . 

It  can  be  shown  using  the  method  of  integral  relations  that  the 
rate  of  flow  Qw  is  given  exactly  by  (Polubarinova- Kochi na  [ 1962,  p.  283], 
Hantush  [1964,  p.  362],  Bear  [1972,  p.  368]): 


- h )/ln(r  /r  ) . 
w e w 


(1) 


Hantush  [1964,  p.  362]  states  that  this  formula  agrees  with  experimental 
results  to  within  2%  . 

The  existence  and  uniqueness  of  flow  towards  a well  has  been 

studied  by  Mauersberger  [1965c] . A variational  principle  has  been  derived 
by  Mauersberger  [1965b]  . 


Numerical  solutions  have  been  obtained  by  many  authors  and 
are  summarized  in  Table  1 below. 


Author  (Method) 

r 

w 

r 

e 

h 

w 

h 

e 

h 

s 

Babbitt  and  Caldwell 

[1948,  p.  27]  (trial- 

.083 

1.0  . 

05875 

. 235 

. 1 128 

free -boundary:  analog) 

and  eleven  other  cases 

Yang  [ 1949]  (trial- 

20 

2,020 

400 

1,600 

1,240  ± 

5 

free -boundary: 

100 

4,  100 

800 

1,600 

1,060  ± 

5 

finite  differences) 

20 

4,020 

1,200 

1,600 

1,320  ± 

10 

20 

4,020 

800 

1,600 

-H 

in 

O' 

10 

20 

4,020 

400 

1,600 

1,  135  ± 

10 

20 

4,020 

0 

1,600 

1,115* 

10 
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Author  (Method) 


w 

Boulton  [ 19  51]  and  50 

J.  A.  Murray  [i960]  200 

(trial-free-boundary;  200 

finite  differences)  1,000 

£ 

6.400 
10,  400 

10.400 
16,000 

W 

3,  200 
3, 200 
0 
0 

6,760 
7,  165 
7.120 
4,000 

s 

5,725 
5,  475 
5,  153 
1, 500 

Kashef,  Touloukhnn,  100 

and  Eadum  [ 1952,  p.73] 

(trial  free  boundary; 
finite  differences).  See 
also  Kashef  [ 1953]. 

4,  100 

800 

1,600 

* 

1,040 

H.  Schmidt  [ 19  56]  1,000 

30,000 

10,000 

8,000 

8,152 

(time-dependent;  1,000 

30,000 

10,000 

6,000 

6,695 

finite  differences)  1,000 

30,000 

10,000 

4,000 

5,71  1 

1,000 

30,000 

10,000 

2,000 

5,  168 

1,000 

30,000 

10,000 

0 

4,894 

H.  P.  Hall  [ 1955]  • 4.8 

76.8 

36 

48 

Graphical 

(trial-free-boundary;  4.8 

76.8 

24 

48 

(see  remarks) 

finite  differences)  4.8 

76.8 

12 

48 

Herbert  [ 1968 , p.  260]  11 

(time-dependent;  resis- 
tance network) 

99 

35 

70 

54.  5 

G.  S.  Taylor  and  Luth in  4.8 
[1969]  (time-dependent; 
finite  differences) 

76.  8 

12 

48 

Graphical 

Neuman  and  Witherspoon  4.8 
[1970]  (trial-free-bound- 
ary;  finite  elements) 

76.8 

12 

48 

Graphical 
(see  remarks) 

Neuman  and  Witherspoon  4.8 
[1971a]  (time-dependent; 
finite  elements) 

76.8 

12 

48 

Graphical 

France,  Parekh,  Peters,  11 
and  Taylor [ 1971,  p.  176 J 
(time-dependent;  finite 
elements) 

99 

35 

70 

* 

46 

Estimated  by  us  from  a graph 

T-'ble  1:  Axisymmetric  flow  towards  a fully  penetrating  well 
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As  can  hr  seen  from  Table  1,  a large  number  of  different  cases 

have  been  considered,  so  that  it  is  difficult  to  compare  the  results.  It  is 

also  unfortunate  that  some  of  the  results  are  given  in  graphical  form. 

However,  the  case  r = 4.8,  r - 76.8,  h = 12,  h = 48  has  been 

considered  by  Hall  [ 19  55]  , G.  S.  Taylor  and  Luthin  [1969]  , and 

Neuman  and  Witherspoon  [ 1970],  and  Neuman  and  Witherspoon  [1970]  give 

a graphical  comparison  of  these  results.  The  results  are 

not  strictly  comparable,  because  Hall  assumed  capillarity  (h  = 3.6) 

c 

and  a lined  well  (h  = 1.0),  Taylor  and  Luthin  assume  partially  saturated 
max 

flow,  and  Neuman  and  Witherspoon  assume  saturated  flow  without  capillarity, 
but  the  results  are  quite  close.  France,  Parekh,  Peters,  and  Taylor  [1971] 
give  a graphical  comparison  of  their  results  with  those  of  Herbert  [1968]. 

R.  L.  Taylor  and  Brown  [1967]  give  graphical  results  which  were 
obtained  using  a general  program.  This  program  is  listed  by  R.  L.  Taylor 
[1966],  and  an  improved  version  is  listed  by  Kealy  and  Busch  [1971]. 
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There  are  several  approximate  analytic  solutions  in  the  literature. 

In  a series  of  papers,  Mauersberger  has  obtained  approximate  solutions 
using  the  variational  method  of  Trefftz  (Mauersberger  [ 1967,  1968,  1968a]  ), 
collocation  (Mauersberger  [ 1968 , 1968b,  1968c]),  least  squares 
(Mauersberger  [1968,  1968c]),  Galerkin's  method  (Mauersberger  [ 1968, 
1968c]),  the  "partition  method"  (Mauersberger  [ 1968  , 1968c]).  Kirkham 
[ 1964]  has  obtained  an  approximate  solution  by  introducing  a "fictitious 
flow  region"  and  using  collocation.  Kirkham  does  not  give  specific 
numbers  and  checks  his  solution  by  comparing  it  with  the  problem  of 
i flow  through  a dam. 

Brutsaert,  Breitenbach,  and  Sunada  [ 1971  ] solve  numerically  the 
corresponding  time-dependent  multiphase  problem. 

3.  4. 1.1. 2.  Fully  penetrating  well:  the  general  case 

The  case  which  has  been  studied  the  most  Intensively  Is  the 
case  of  N parallel  aquifers  (Figure  1):  aquifer  i is  of  depth  a.  and 
has  permeability  k^.  (Note  that  the  ordering  of  the  aquifers  is  sometimes 
reversed,  aquifer  1 being  at  the  bottom.  ) In  other  words, 

k = k(y)  = k , if  Z a <y<  a . (1) 

j=i-l  J j=i  ] 
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Neumann  and  Witherspoon  ( 1969]  give  a detailed  discussion  of 
transient  flow  for  a multiple  aquifer  (Hyure  1)  under  the  assumption  that  the 
FH  is  horizontal. 

Mauersberger  [ 1969]  shows  that  it  is  possible  to  use  the  method 
of  integral  relations  to  obtain  an  expression  for  the  rate  of  flow  Qw  when 
the  permeability  k is  of  the  form 

k = k(x,  0,  y)  = kj(x)k2(0,  y) 

where  (x,  0,  y)  are  axial  cylindrical  coordinates.  This  includes  the  cases 
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k = k (x)  and  k = k (y)  as  special  cases.  In 'particular,  Mauersberger 
1 2 

qives  the  value  of  Q when  k is  of  the  form  (1). 

w 

Youngs  [ 1971b I considers  the  case  when  k = k(y)  and  when 
precipitation  occurs  at  a rate  R per  unit  area.  For  the  special  case 


r 


k(y)  = < 


k2,  0 < z < a^, 


*V  a2  - Z’ 


h = a + a , and  h = a.  +1  + H',  Youngs  obtains  the  formula 

w 1 2 e 1 2 

Qw  (H')2k1R  2RH'(kj3j  + kza2)  R(  <re/rw>  -1) 

’ Rr 2 t n{r  /r  ) + RrLnfr/r  ) 

w w e w w e w 

If  there  is  no  precipitation  (R  = 0)  then  this  formula  of  Youngs  is  of 


course  contained  in  the  results  of  Mauersberger. 


Lewis  and  Uuinphcson  [ 1 ‘>7 3 ] use  finite  elements  to  solve  the  time- 
dcpendenL  problem  of  a well  subjected  to  electrical  forces  (figure  2). 
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Figure  2. 


Free  surface  profiles  under  varying  voltage 
gradients  after  2 1/2  davs  (based  on  Lewis 
and  Humpheson  [197  3,  p.  612]) 


The  geometry  is  shown  in  Figure  1. 


Numerical  solutions  have  been  obtained  by  several  authors:  Boreli 
1 1955]  (trial-free-boundary;  finite  differences);  R.  L.  Taylor  and  Brown 
[1967]  using  the  program  of  R.  L.  Taylor  [ 1966]  (trial-free-boundary; 
finite  elements);  Cooley  f 1 97 1 ] (time-dependent  partially-saturated  flow; 
finite  differences).  Unfortunately,  only  Boreli  states  the  dimensions 
clearly,  and  all  the  authors  give  their  results  in  graphical  form. 


Breitenodcr  [ 1942,  p.  84  and  p.  86]  considers  the  plane  problem 

for  d = r = 00  . 
e 

It  would  seem  that  the  problem  deserves  further  study. 

3.  4. 1.3.  Multiple  wells 

The  problem  of  multiple  wells  is  a three-dimensional  problem 

on  which  as  yet  little  work  has  been  done. 

Youngs  [ 1971b]  considers  the  approximate  problem  of  a single 

well  in  a multiple  aquifer  (see  Figure  1 of  section  3.  4. 1.1.  2)  for  which 

there  is  no  flow  at  r = r but  precipitation  R per  unit  area.  Using 

e 

integral  identities  Youngs  obtains  upper  and  lower  bounds  for  the  drawdown 


for  a three-layered  soil. 


4 . Two -fluid  tlows 


Two-fluid  porous  flow  TBS  arise  when  a lighter  fluid  (fluid  1)  lies 
on  top  of  a heavier  fluid  (fluid  Z ) , the  interface  between  the  fluids  berng 
a FB.  In  the  analysis  of  such  problems  it  is  usually  assumed  that  the 
lower  fluid  is  at  rest.  The  boundary  conditions  on  the  interface  are 
given  in  section  2.2. 

Two-fluid  porous  flow  FBPS  arise  in  two  important  contexts  namely 
sea  water/fresh  water  FBPS  in  coastal  areas  and  water/gas/oil  FBPS  in 
the  oil  industry,  and  these  are  discussed  in  the  first  two  subsections  below. 
The  final  subsection  discusses  the  FBPS  which'arise  when  one  considers  the 
FB  at  the  microscopic  level. 

4.1.  Salt  water  - fresh  water  interfaces 

In  coastal  regions  there  is  an  interaction  between  fresh  water  being 
applied  to  the  ground  surface  by  rainfall  or  irrigation,  and  salt  water  from 
the  sea.  The  fresh  water  is  lighter  and  forms  a layer  on  top  of  the  salt  water. 

It  is  of  importance  to  prevent  the  encroachment  of  the  sea  water  because 
the  salt  destroys  agricultural  land  and  contaminates  water  supplies. 

In  studying  salt  water/fresh  water  problems  it  is  usual  to  assume  that 
the  salt  water  is  at  rest,  and  this  assumption  will  be  made  throughout  unless 
explicitly  stated  otherwise.  Indeed,  sc  far  as  we  are  aware  only  Childs  [ 1950] 
has  considered  the  case  when  the  salt  water  is  in  motion.  In  this  connection 
it  is  of  interest  that  the  observed  fresh  water/salt  water  interface  in  the  Miami 
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area  shows  that  motion  of  the  sea  water  must  occur  (Cooper,  Kohout,  Henry,  and 
Glover  [ 1964,  p.  C12]  ). 

Of  course,  the  assumption  of  a fresh  water/salt  water  interface  is  only 
an  approximation,  and  in  the  general  case  it  is  necessary  to  take  account  of 
the  diffusion  of  the  salt. 


4.1.1.  Salt  water  - fresh  water:  the  Ghyben-Herzberq  lens 


The  first  salt  water  - fresh  water  FBP  to  be  studied  occurs  when 


fresh  water  is  provided  by  either  rainfall  or  irrigation  to  an  island  (axisym- 
metric)  or  isthmus  (plane).  A bubble  of  fresh  water,  often  called  the 
Ghyben-Horzbcrg  lens,  forms  underneath  the  island  (figure  1).  There 


are  two  FBS:  the  fresh  water/air  interface  and  the  fresh  water/salt  water 
interface. 


Childs  [ 1950]  used  the  trial-free-boundary  method  with  conducting 
paper  to  obtain  the  solution  for  the  case  of  a rectangular  isthmus. 

Childs  allows  the  salt  water  to  be  at  different  levels,  which  can  arise 
because  of  tidal  effects  (Figure  2).  Youngs  [1971,  1971a]  obtains  bounds 
for  the  rate  of  flow  with  and  without  a well  on  the  island. 
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Figure  2.  The  Ghyben-Herzberg  lens  for  a 

rectangular  Isthmus  (based  on  Childs 
f 1950.  p.  177]  ) 


Ackerman  and  Intong  [ 1966]  consider  a related  problem  in  which 
fresh  water  is  supplied  from  a canal  with  impervious  walls.  It  is  assumed 
that  the  upper  surface  is  known  and  that  the  salt  water  is  at  rest  (Figure  3). 
An  analytical  solution  is  found  using  the  hodograph  method. 
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Figure  3.  Fresh  water  lens  fed  from  a canal 


Charmonman  [ 1966]  considers  the  problem  shown  in  Figure  4: 
fresh  water  is  supplied  to  a canal  in  the  center  of  an  isthmus  with  an 
impervious  ground  surface.  It  is  assumed  that  the  salt  water  is  at 
rest.  This  problem  was  solved  by  Charmonman  [ 1966]  using  finite 
differences  in  the  4>4>-plane  and  by  Mason  and  Farkas  [ 1971,  1972]  using 
the  trial-free-boundary  method  with  superposition. 

Keuning  [ 1967]  considers  the  abstraction  of  fresh  water  from  an 
underground  well  (sink)  in  a nonhomogeneous  circular  island  fed  by  rain 
water.  Keuning  makes  the  Dupuit  assumption  that  the  flow  is  parallel  to 
the  surface,  and  it  would  be  of  interest  to  consider  the  exact  problem. 
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Figure  4.  Canal  with  impervious  surface 


4.1.2.  Coastal  aquifers 

The  geometry  is  shown  in  Figure  1:  fresh  water  flows  from  a coastal 
aquifer  into  the  sea.  There  are  two  FBS  namely  the  upper  fresh  water/air  FB 
and  the  lower  fresh  water/snlt  water  FB.  Charmonman  [ 1965]  uses  the 
hodograph  method  to  obtain  an  analytical  solution  for  the  case  when  the 
outflow  surface  BC  is  horizontal. 

Henry  [1959]  uses  the  hodograph  method  to  obtain  the  exact  solutions 
for  flow  from  a horizontal  aquifer  of  finite  depth  with  vertical  outflow  face 
(Figure  2a)  and  horizontal  outflow  face  (Figure  2b).  (The  problem  of 
Figure  2b  for  the  case  of  infinite  depth  was  solved  earlier  by  Glover  [ 1959]  ). 
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Figure  1.  Flow  from  a coastal  aquifer 


The  coupled  equations  for  Darcy  flow  of  water  and  dispersion  of  salt 
have  been  solved  numerically  for  the  problem  of  Figure  2a  by  Ilenry  (see 
Cooper,  Kohout,  Henry,  and  Glover  [ 1964,  p.  C70]  ) and  Pinder  and 
Cooper  [ 1970]  . 
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(a)  vertical  outflow 


(b)  horizontal  outflow 


Figure  2.  Horizontal  aquifers  of  finite  depth 


Dejong  [ 1965]  uses  the  hodograph  method  to  obtain  the  analytical 
• Int io’  for  flow  in  a semi-infinite  aquifer  with  a surface  sink  (Figi 3). 


i)n  Jong  alr.o  gives  experimental  results  which  are  in  e \ < '.!•  . * uqn  ent 
with  the  theoretical  results. 
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Figure  3.  Semi-infinite  aquifer  with  sink 

Baiocchi,  Comincioli,  Magenes,  and  Pozzi  [ 1973,  p.  8 and  p.  39] 

consider  the  generalization  of  the  problem  of  Figure  2a  obtained  by 

allowing  a fresh  water/air  interface  (Figure  4).  It  is  assumed  that  the 

aquifer  is  so  wide  that  FG  is  a line  of  equipotcntial.  Baiocchi, 

Comincioli,  Magenes,  and  Pozzi  reformulate  the  problem  as  a variational 

inequality  and  prove  existence  and  uniqueness  of  the  solution.  Baiocchi, 

Comincioli,  Guerri,  and  Volpi  [ 1973,  p.  29]  use  variational  inequalities 

with  finite  differences  to  obtain  a numerical  solution.  Youngs  [1971a]  derives 
the  exact  value  of  the  rate  of  flow. 

In  some  coastal  aquifers  the  aquifer  is  divided  into  several  strata  by 
interlying  semi-pervious  layers.  As  a result,  seawater  intrudes  into  each 
of  the  separate  aquifers  and  forms  a series  of  wedges.  Such  coastal 
aquifers  occur  near  New  York,  and  off  the  coasts  of  Israel  and  the 
Netherlands  (M.  A.  Collins  and  Gelhar  [ 1971  ] ) . The  geometry  is 


shown  in  Figure  5.  Collins  and  Gelhar  [1971]  give  an  approximate 
analysis  based  upon  the  Dupuit  approximation. 


Figure  5.  A layered  aquifer 


4.1.3.  Land  reclamation 


In  certain  areas  of  the  world,  such  as  the  central  plain  of  Thailand, 
it  is  conceivable  that  it  will  be  possible  to  reclaim  land  from  the  sea  by 
applying  fresh  water  to  the  surface. 

Charmonmon  [ 1967]  uses  finite  differences  in  the  <}>4>-plnne  to 
analyse  an  array  of  parallel  drains  and  canals  (Figure  la).  Ackerman  and 
Chang  [ 1971]  use  the  hodograph  method  to  obtain  an  analytical  solution  for 
the  case  when  fresh  water  is  applied  to  the  surface  and  then  (to  reduce  the 
loss  of  fresh  water)  pumped  from  a series  of  drains  (Figure  lb). 
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Figure  1.  Methods  of  land  reclamation 
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4.2.  Other  configurations 


4.2.1.  Up-conlnq 

When  a layer  of  one  fluid  (fluid  1)  lies  over  n layer  of  a second 
fluid  (fluid  2)  and  fluid  1 is  pumped  out  from  one  or  more  wells  (or  sinks) 
then  the  interface  between  the  fluids  moves  towards  the  wells  and  this 
is  called  up-coninq  because  the  interface  often  approximately  assumes 
the  shape  of  a cone.  If  the  rate  of  pumping  is  too  high  then  fluid  2 may 
enter  the  well  and  contaminate  the  well  output.  Three  especially 
important  cases  occur  in  practice: 


Fluid  1 

Fluid  2 

Special  terminology 

oil 

water 

water-coning 

gas 

oil 

gas-coning 

fresh  water 

salt  water 

Water-coning  is  discussed  by  Muskat  [ 19  37  , p.  480,  Muskat  [1949], 
Pirson  [ 19  58] ; gas -coning  is  discussed  by  Muskat  [ 19  37 , p.  689]: 
and  up-coning  is  discussed  by  Bear  [197  2,  p.  553  and  p.  569]. 

The  hodograph  method  has  been  used  to  obtain  exact  solutions 
for  two  cases:  (i)  a single  sink  when  fluids  1 and  2 are  both  semi- 
infinite (Figure  1 and  Bear  [1972,  p.  553]);  (ii)  a line  of  sinks  when 
fluid  1 is  semi-infinite,  fluid  2 is  of  finite  depth,  and  the  flow  is  critical, 
that  is,  the  FB  has  a cusp  (Figure  2 and  Kidder  [19  56]). 
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Figure  1:  A single  sink  and  two  semi-infinite  fluids 
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Figure  2:  Critical  flow  for  n line  of  sinks  with  fluid  1 
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Karplur,  | 19  56]  describes  the  use  of  the  trial-freo-boundary 
method  with  resistance  networks  for  a single  axisymmetric  well  of  finite 
dimensions  (Figure  3).  Karplus  [1956,  p.  244]  states  that  he  has 
investigated  several  problems  but  does  not  give  his  results. 


Figure  3:  Water  cone  in  vicinity  of  oil  well 
(based  on  Karplus  [1956,  p.  241]) 

Finally,  we  observe  that  Welge  and  Weber  [1964]  solve 
numerically  the  axisymmetric  time -dependent  equations  for  two  immiscible 
fluids  and  obtain  good  agreement  between  their  results  and  experimental 
and  field  observations. 

*.2.2.  Inclined  reservoirs 

If  oil  is  present  in  an  inclined  stratum  then  FBPS  of  the  form 
shown  in  Figure  1 arise  (Bear  (1972,  p.  534]  ).  Superficially  at  least  such 
problems  bear  a similarity  to  fluid  mechanics  FBPS  for  infinite  bubbles. 


Figure  1 . An  inclined  plane  reservoir 
4.3.  Interfacial  instability 

When  a viscous  fluid  1 in  a porous  medium  is  driven  forwards 

by  the  pressure  of  another  driving  fluid  2,  the  interface  between  the 

two  fluids  is  liable  to  be  unstable  if  p < p.  where  p.,  p are  the 

1 1 z 

viscosities  of  the  two  fluids.  Heuristicully , one  can  see  why  such 
instability  occurs:  the  driving  fluid  (fluid  2)  encounters  less  resistance 
than  fluid  1 so  that  for  a given  driving  pressure  the  flow  rate  will  be 
greater  if  the  interface  is  distorted  so  that  fluid  2 can  flow  leaving  some 
of  fluid  1 behind. 

Saffman  and  Taylor  [ 19  58]  modeled  the  problem  using  the  analogy 
between  porous  flow  and  viscous  flow  in  a Hele-Shaw  cell  (see 
Bear  [ 1972,  p.  687  |).  Figure  1 shows  how  an  interface,  which  was 
originally  straight,  deforms  into  a number  of  fingers,  and  for  this  reason 
the  phenomen  is  called  fingering- 
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The  phenomenon  o[  fingering  con  occur  in  oil  production  in  which 
a more  viscous  fluid  (oil)  is  being  driven  by  a less  viscous  fluid  (water). 
It  seems  to  us,  however,  that  there  is  some  confusion  in  the  literature 
between  fingering  and  water-coning.  Figure  2 shows  a line  of  oil  wells 
producing  oil. 


tS  1/ 


The  oil/water  interface  is  distorted  in  the  neighborhood  of  each  well, 

<.;-d  ■ are  formed.  However,  the  mechanism  of  the  phenomena 

. ■ -lines  1 and  2 is  quite  different:  in  Figure  1 the  driving 
pressure  is  uniform  and  the  fingers  are  caused  by  the  instability  of  the 
surface;  in  Figure  2 the  "fingers"  correspond  to  the  spacing  of  the  wells. 

It  thus  seems  to  us  that  it  is  incorrect  to  compare  fingering  and  water- 
coning  as  done  by  Saffman  and  Taylor  [ 19  58,  p.  314]  . Water-coning  is 
discussed  in  section  4.2. 

To  analyse  fingering,  Saffman  and  Taylor  [ 19  58 , p.  318] 
consider  the  case  when  there  is  an  infinite  set  of  equal  and  equally 
spaced  fingers  all  advancing  at  the  same  speed.  Since  each  finger  is 
then  identical  mathematically  with  all  the  others  and  the  fluid  on  the 
straight  lines  halfway  between  neighbors  has  no  transverse  component 
of  velocity,  it  suffices  to  consider  a single  finger  propagating  itself 
in  a channel  of  fixed  width. 

Saffman  and  Taylor  use  the  hodograph  method  to  obtain  an  analytical 
solution,  or  rather  a family  of  solutions,  because  the  solution  contains 
the  parameter 

^ width  of  channel 

asymptotic  width  of  finger 

Experimentally,  Saffman  and  Taylor  find  that  X is  very  close  to 
1/2,  but  they  are  unable  to  explain  this:  the  maximum  velocity  occurs 
when  X = 0;  the  maximum  and  minimum  rate  of  energy  dissipation  occur 
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when  X.  = 0 or  X.  = 1;  a stability  analysis  shows  that  the  motion  is 
unstable  for  all  X.. 

Saffman  [1959]  obtains  a family  of  exact  solutions  of  the  time- 
dependent  problem  which  begin  as  an  approximately  sinusoidal  perturbation 
and  grow  into  fingers.  Saffman  [ 1959,  p.  159]  observes  that  the  growth 
is  slowest  when  X.  =r  1/2  but  adds  that  the  physical  significance  of  this 
is  not  clear.  Similar  exact  solutions  were  obtained  by  Jacquard  and 
Seguier  [ 1962]  . 

Wo  conclude  with  some  observations: 

(i)  Bear  [1972,  p.  514]  discussed  fingering  and  gives  references. 

The  approach  followed  is  different  from  that  of  Saffman  and  Taylor, 
and  emphasizes  the  conditions  for  stability  or  instability  rather  than 
the  shape  of  the  fingers. 

(ii)  There  is  considerable  similarity  between  fingering  and  other 
interfacial  instability  phenomena  such  as: 

(a)  Dendrite  formation  in  crystals. 

(b)  Lubrication  cavitation  in  bearings . 

(c)  Taylor  instability  at  the  accelerated  interface  between 
fluids  of  different  densities. 

(iii)  As  Saffman  and  Taylor  observe,  the  non-uniqueness  of  fingers 
(because  of  the  arbitrary  parameter  X.)  is  similar  to  the  non-uniqueness 
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of  bubbles  moving  in  tubes.  In  both  cases  it  would  seem  that  an 
additional  physical  condition  must  be  imposed. 

(iv)  S.  Richardson  [ 1972]  considers  similar  problems  arising  in  the 
injection  molding  of  plastics. 

(v)  The  analysis  of  Saffman  and  Taylor  is  for  plane  problems;  so 
far  as  we  are  aware,  axisymmetric  problems  have  not  been  considered. 


S.  Coupled -field  problems 


Under  coupled-field  porous  flow  FBPS  we  understand  porous  flow 
FBPS  which  involve,  in  a non-trivial  way,  one  of  the  other  fields  of 
continuum  mechanics  such  as  heat  flow  on  electromagnetism. 

There  are  many  possible  coupled-field  FBPS  involving  dispersion 
or  diffusion.  Bear  [ 1972,  chapter  10]  gives  an  excellent  discussion  and 
lists  the  following  applications:  (a)  the  transition  zone  between  saltwater 
and  fresh  water  in  coastal  aquifers;  (b)  artificial  recharge  operations  where 
water  of  one  quality  is  introduced  into  aquifers  containing  water  of  a 
different  quality  (mixing  of  water  due  to  hydrodynamic  dispersion  is  among 
the  various  objectives  of  artificial  replenishment);  (c)  secondary  recovery 
techniques  in  oil  reservoirs,  where  an  injected  fluid  dissolves  the 
reservoir's  oil;  (d)  radioactive  and  reclaimed  sewage  waste  disposal  into 
aquifers;  (e)  the  use  of  tracers,  such  as  dyes,  electrolytes  and  radioactive 
isotopes,  in  hydrology,  petroleum  engineering  and  other  scientific  and 
engineering  research  projects;  (f)  the  use  of  reactors  packed  with  granular 
material  in  the  chemical  industry;  and  (g)  the  movement  of  fertilizers  in 
the  soil  and  the  leaching  of  salts  from  the  soil  in  agriculture.  We  have 
not  searched  the  literature  for  FBPS  involving  these  phenomena,  although 
such  problems  doubtless  occur. 

Another  important  coupled-field  problem  involves  heat  transfer. 
Kirkham  and  Powers  [ 197  2,  p.  462]  observe  that  soil  temperature  affects 
the  germination  of  seeds,  the  chemical  reactions  which  release  nutrients, 


and  the  availability  of  water.  Also  (Kirkham  and  Powers  [ 1972,  p.  483]) 
the  freezing  ol  soil  is  of  importance  in  agriculture  and  in  civil  engineering 
projects  involving  the  building  of  roads  and  pipelines.  Bear  [ 1972,  p.  641  ] 
discusses  heat  transfer  in  a porous  medium.  Here  again,  we  have  not 
searched  the  literature  for  FBPS. 

A third  coupled-field  problem  arises  in  connection  with  electro- 
kinetic  effects:  the  movement  of  water  under  the  influence  of  an  electric 
field.  This  is  duscussed  further  in  section  1 . 6 and  the  flow  towards  a 
well  has  been  treated  by  Lewis  and  Humpheson  [ 1973]  (see  section  3.  4.1.1.  2). 
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